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Abstract 

Polyurea  is  a  special  type  of  elastomer  that  features  fast  setting  time  as  well  as  good  chemical 
and  fire  resistance.  It  has  also  good  mechanical  properties  such  as  its  high  toughness-to-density 
ratio  and  high  strain  rate-sensitivity,  so  its  application  is  recently  extended  to  structural  purpose 
to  form  sandwich- type  or  multi-layered  plates.  Those  structures  can  be  used  for  retrofitting  of 
military  vehicles  and  historic  buildings,  absorbing  energy  during  structural  crash. 

In  order  to  investigate  its  behavior  of  hysteresis  as  well  as  rate-sensitivity,  three  different 
testing  systems  are  used  to  cover  a  wide  range  of  strain  rates  up  to  strain  of  100%.  In  view 
of  impact  and  blast  events,  the  virgin  state  of  polyurea  is  considered  throughout  the  exper¬ 
iments.  First,  a  hydraulic  universal  testing  machine  is  used  to  perform  uniaxial  compressive 
loading/unloading  tests  in  order  to  investigate  its  hysteresis  behavior  at  low  strain  rates  (0.001/ s 
to  10/s).  Second,  two  distinct  gas-gun  split  Hopkinson  pressure  bar  [SHPB]  systems  are  em¬ 
ployed  to  cover  high  strain  rates:  a  nylon  bar  system  (700/s  to  1200/s)  and  an  aluminum  bar 
system  (2300/s  to  3700/s).  Lastly,  the  rate-sensitivity  for  intermediate  strain  rates  (10/s  to 
1000/s)  is  characterized  using  a  modified  SHPB  system.  The  device  is  composed  of  a  hydraulic 
piston  along  with  nylon  input  and  output  bars. 

A  finite  strain  constitutive  model  of  polyurea  is  presented  in  order  to  predict  the  hysteresis 
and  rate-sensitivity  behavior.  The  1-D  rheological  concept  of  two  Maxwell  elements  in  parallel 
is  employed  within  the  framework  of  the  multiplicative  decomposition  of  the  deformation  gra¬ 
dient.  Model  parameters  are  calibrated  based  on  the  uniaxial  compressive  tests  at  various  rates. 
The  corresponding  algorithms  is  implemented  as  a  user-defined  material  subroutine  VUMAT 
for  ABAQUS/Explicit,  and  used  to  predict  the  response  of  polyurea.  The  proposed  consti¬ 
tutive  model  reasonably  captures  the  experimentally  observed  asymmetric  rate-sensitivity  and 
stress-relaxation  behavior:  strong  rate-sensitivity  and  large  amount  of  stress  relaxation  dur¬ 
ing  loading  phase,  but  weak  rate-sensitivity  and  smaller  amount  of  stress  relaxation  during 
unloading  phase.  In  order  to  validate  the  proposed  model,  various  dynamic  punching  tests 
are  performed,  and  their  results  are  well  compared  with  the  model  predictions  during  loading 
although  the  prediction  of  unloading  behavior  can  be  further  improved. 
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Chapter  1 


Introduction 


1.1  Motivation 

Polyurea  has  been  generally  used  for  coating  purpose  due  to  its  good  chemical  properties  such  as 
fast  setting  time  (few  minutes  or  less)  and  good  water/chemical/fire  resistance.  Polyurea  is  used 
on  metallic  substrates  where  it  provides  corrosion  and  abrasion  resistance  in  harsh  environments. 
Moreover,  it  has  also  good  mechanical  properties  such  as  high  toughness  for  its  low  density  and 
high  rate-sensitivity.  Thus,  its  application  is  recently  extended  to  structural  purpose  to  form 
sandwich-type  or  multi-layered  plates,  which  can  be  used  for  retrofitting  of  military  vehicles, 
historic  buildings,  gas/oil  pipelines  and  marine  structures,  absorbing  energy  during  structural 
crash  and  holding  metal/brick  fragments  even  after  structural  failure.  Like  other  elastomers, 
polyurea  is  viscous  material  so  its  mechanical  properties  are  rate  dependent.  Especially  under 
extreme  loading  conditions  such  as  blast,  projectile  and  explosive  loadings,  polyurea  becomes 
very  attractive;  its  lightness  is  good  for  operational  purpose,  and  its  mechanical  strength  is 
enhanced  under  extreme  loading  conditions.  This  structural  purpose  of  polyurea  under  extreme 
loading  conditions  motivates  our  research  on  a  wide  range  of  strain  rates  up  to  large  strains. 
In  view  of  the  simulation  of  blast  and  impact  events,  we  limit  our  attention  to  the  mechanical 
behavior  of  polyurea  in  its  virgin  state. 

Limited  experimental  study  of  virgin  polyurea  has  been  published  on  the  rate-sensitivity 
behavior  for  a  wide  range  of  strain  rates.  Polyurea  shows  a  highly  nonlinear  viscoelastic 
behavior  at  finite  strains  (e.g.  Amirkhizi  et  al.,  2006,  Bogoslovov  and  Roland,  2007,  Roland 
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et  al.,  2007).  The  stress-strain  response  of  most  polymeric  materials  shows  a  pronounced 
strain  rate  sensitivity  at  low,  intermediate  and  high  strain  rates.  Various  authors  published 
experimental  results  on  the  strain  rate  sensitive  response  of  amorphous  glassy  polymers  (e.g. 
Chou  et  al.,  1973,  Boyce  et  al.,  1988,  Walley  et  al.,  1989,  Cady  et  al.,  2003,  Siviour  et  al.,  2005, 
Mulliken  and  Boyce,  2006,  Mulliken  et  al.,  2006),  crystalline  glassy  polymers  (e.g.  Chou  et  al., 

1973,  Bordonaro  and  Krempl,  1992,  Cady  et  al.,2003,  Siviour  et  al.,  2005,  Khan  and  Farrokh, 
2006)  and  elastomers  (e.g.  Gray  et  al.,  1997,  Rao  et  al.,1997,  Song  and  Chen,  2003,  2004,  Hoo 
Fatt  and  Bekar,  2004,  Shergold  et  al.,  2006,  Roland,  2006)  including  polyurea  (e.g.  Amirkhizi 
et  al.,  2006,  Roland  et  al.,  2007,  Sarva  et  al.,  2006).  However,  only  few  experimental  studies 
deal  with  the  intermediate  strain  rate  behavior  of  elastomers  at  large  deformations.  Sarva 
et  al.  (2006)  performed  intermediate  strain  rate  compression  tests  on  polyurea  for  maximum 
strains  greater  than  100%  but  the  strain  rates  were  only  14  ~  80/s.  The  same  research  group 
also  obtained  test  results  for  a  strain  rate  of  800/s  using  a  very  long  aluminum  SHPB  system. 
Roland  et  al.  (2007)  characterized  the  tensile  behavior  of  polyurea  over  a  strain  rate  range  of 
14  ~  573/s  and  up  to  strains  of  more  than  300%. 

As  for  the  experimental  study  on  the  rate-sensitivity  of  polyurea,  little  research  has  been 
reported  on  modeling  of  virgin  polyurea  under  loading/unloading  conditions  for  a  wide  range 
of  strain  rates.  Finite  viscoelasticity  models  of  elastomers  may  be  formulated  using  the  so- 
called  hereditary  integral  approach  (Coleman  and  Noll,  1961,  Bernstein  et  al.,  1963,  Lianis, 
1963,  McGuirt  and  Lianis,  1970,  Leonov,  1976,  Johnson  et  al.,  1994,  Haupt  and  Lion,  2002, 
Amirkhizi  et  al.,  2006)  but  their  validity  is  often  limited  to  a  narrow  range  of  strain  rates 
(Yang  et  al.,  2000,  Shim  et  al.,  2004,  Hoo  Fatt  and  Ouyang,  2007).  As  an  alternative  to  the 
hereditary  integral  approach,  the  framework  of  multiplicative  decomposition  of  the  deformation 
gradient  (Kroner,  1960  and  Lee,  1969)  is  frequently  used  in  finite  viscoelasticity  (e.g.  Sidoroff, 

1974,  Lubliner,  1985,  Le  Tallec  et  al.,  1993,  Reese  and  Govindjee,  1998,  Huber  and  Tsakmakis, 
2000).  In  that  framework,  the  nonlinear  viscoelasticity  of  elastomers  is  commonly  described 
through  a  rheological  spring-dashpot  models  of  the  Zener  type  (e.g.  Roland,  1989,  Johnson 
et  al.,  1995,  Bergstrom  and  Boyce,  1998,  Huber  and  Tsakmakis,  2000,  Quintavalla  and  John¬ 
son,  2004,  Bergstrom  and  Hilbert,  2005,  Qi  and  Boyce,  2005,  Areias  and  Matous,  2008,  Hoo 
Fatt  and  Ouyang,  2008,  Tomita  et  al.,  2008).  As  for  the  hereditary  integral  approach  based 
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models,  however,  most  multiplicative  decomposition  based  models  have  been  also  experimen¬ 
tally  validated  for  a  narrow  range  of  strain  rates.  Quintavalla  and  Johnson  (2004)  adopted  the 
Bergstrom-Boyce  model  to  describe  the  dynamic  behavior  of  cis-(l,4)  polybutadiene  at  high 
strain-rate  of  3000/s  to  5000/s.  Recently,  Hoo  Fatt  and  Ouyang  (2008)  proposed  a  thermo¬ 
dynamically  consistent  constitutive  model  with  modified  neo-Hookean  rubber  elastic  springs 
to  describe  the  steep  initial  stiffness  of  virgin  butadiene  rubber  under  tensile  and  compressive 
loading  at  intermediate  strain  rates  (76/s  to  450/s). 

1.2  Objective  and  Tasks 

The  objective  of  the  present  dissertation  is  to  develop  a  finite  constitutive  model  of  virgin 
polyurea  for  a  wide  range  of  strain  rates  under  loading/unloading  condition.  The  material  is 
assumed  to  be  isotropic,  and  the  isothermal  condition  is  considered  at  the  room  temperature. 
In  order  to  achieve  the  objective,  the  following  tasks  are  set: 

•  To  characterize  the  mechanical  properties  of  virgin  polyurea  at  low,  intermediate  and  high 
strain  rates  for  large  strains; 

•  To  develop  a  finite  constitutive  model  of  virgin  polyurea; 

•  To  demonstrate  the  validity  of  the  proposed  constitutive  model  by  performing  dynamic 
punching  tests  and  simulations. 

1.3  Outline  of  Dissertation 

This  dissertation  is  composed  of  five  chapters.  Excluding  Chapter  1  and  Chapter  5,  each 
chapter  addresses  one  specific  topic  and  it  is  self-contained  because  it  has  already  been  published 
or  submitted  for  publications.  A  list  of  publications  related  to  this  dissertation  is  presented  in 
Appendix  B. 

Chapter  2  presents  experimental  results  on  polyurea.  The  strain  rate  sensitivity  of  polyurea 
is  characterized  using  a  modified  split  Hopkinson  pressure  bar  (SHPB)  system.  The  device  is 
composed  of  a  hydraulic  piston  along  with  nylon  input  and  output  bars.  In  combination  with 
an  advanced  wave  deconvolution  method,  the  modified  SHPB  system  provides  an  unlimited 


14 


measurement  time,  and  thus  can  be  used  to  perform  experiments  at  low,  intermediate  and  high 
strain  rates.  A  series  of  compression  tests  of  polyurea  is  performed  using  the  modified  SHPB 
system.  In  addition,  conventional  SHPB  systems  as  well  as  a  universal  hydraulic  testing  machine 
are  employed  to  confirm  the  validity  of  the  modified  SHPB  technique  at  low  and  high  strain 
rates.  The  analysis  of  the  data  at  intermediate  strain  rates  shows  that  the  strain  rate  is  not 
constant  due  to  multiple  wave  reflections  within  the  input  and  output  bars.  It  is  demonstrated 
that  intermediate  strain  rate  SHPB  experiments  require  either  very  long  bars  (>  20m)  or  very 
short  bars  (<  0.5m  )  in  order  to  achieve  an  approximately  constant  strain  rate  throughout  the 
entire  experiment. 

Chapter  3  is  devoted  to  a  constitutive  model  for  polyurea  based  on  the  experimental  results. 
Continuous  loading  and  unloading  experiments  are  performed  at  different  strain  rates  to  char¬ 
acterize  the  large  deformation  behavior  of  polyurea  under  compressive  loading.  In  addition, 
uniaxial  compression  tests  are  carried  out  with  stair-like  strain  history  profiles.  The  analysis 
of  the  experimental  data  shows  that  the  concept  of  equilibrium  path  may  not  be  applied  to 
polyurea.  This  finding  implies  that  viscoelastic  constitutive  models  of  the  Zener  type  are  no 
suitable  for  the  modeling  of  the  rate  dependent  behavior  of  polyurea.  A  new  constitutive  model 
is  developed  based  on  a  rheological  model  composed  of  two  Maxwell  elements.  The  soft  rub¬ 
bery  response  is  represented  by  a  Gent  spring  while  nonlinear  viscous  evolution  equations  are 
proposed  to  describe  the  time-dependent  material  response.  The  eight  material  model  parame¬ 
ters  are  identified  for  polyurea  and  used  to  predict  the  experimentally-measured  stress-strain 
curves  for  various  loading  and  unloading  histories.  The  model  provides  a  good  prediction  of 
the  response  under  monotonic  loading  over  wide  range  of  strain  rates,  while  it  overestimates 
the  stiffness  during  unloading.  Furthermore,  the  model  predictions  of  the  material  relaxation 
and  viscous  dissipation  during  a  loading/unloading  cycle  agree  well  with  the  experiments. 

Chapter  4  presents  the  validation  application  for  the  proposed  constitutive  model.  Punch 
indentation  experiments  are  performed  on  10mm  thick  polyurea  layers  on  a  steel  substrate. 
A  total  of  six  different  combinations  of  punch  velocity,  punch  size  and  the  lateral  constraint 
conditions  are  considered.  Furthermore,  the  time  integration  scheme  for  a  newly-developed  rate- 
dependent  constitutive  material  model  is  presented  and  used  to  predict  the  force-displacement 
response  for  all  experimental  loading  conditions.  The  comparison  of  the  simulations  and  the 
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experimental  results  reveals  that  the  model  is  capable  to  predict  the  loading  behavior  with  good 
accuracy  for  all  experiments  which  is  seen  as  a  partial  validation  of  the  model  assumptions 
regarding  the  pressure  and  rate  sensitivity.  As  far  as  the  unloading  behavior  is  concerned,  the 
model  predicts  the  characteristic  stiff  and  soft  phases  of  unloading.  However,  the  comparison 
of  simulations  and  experiments  also  indicates  that  the  overall  model  response  is  too  stiff.  The 
results  from  cyclic  compression  experiments  suggest  that  the  pronounced  Mullins  effect  needs 
to  be  taken  into  account  in  future  models  for  polyurea  to  improve  the  quantitative  predictions 
during  unloading. 

Chapter  5  summarizes  the  main  contributions  of  the  dissertation,  and  presents  suggestions 
for  future  studies. 
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Chapter  2 


Experimental  Work 


2.1  Introduction 

Polyurea  is  a  special  type  of  elastomer  which  is  widely  used  as  coating  material.  It  features  a 
fast  setting  time  (few  minutes  or  less)  as  well  as  good  chemical  and  fire  resistance.  Polyurea 
is  frequently  used  on  metallic  substrates  where  it  provides  corrosion  and  abrasion  resistance 
in  harsh  environments.  Applications  include  transportation  vehicles,  pipelines,  steel  buildings 
or  marine  constructions.  More  recently,  polyurea  is  also  considered  for  the  blast  protection 
of  transportation  vehicles  because  of  its  high  toughness-to-density  ratio,  in  particular  at  high 
strain  rates.  It  is  the  objective  of  this  work  to  characterize  the  mechanical  properties  of  polyurea 
at  low,  intermediate  and  high  strain  rates. 

The  mechanical  properties  of  most  metallic  engineering  materials  exhibit  only  a  weak  rate- 
dependence  at  strain  rates  below  100/s.  Therefore,  metals  are  usually  tested  either  at  very 
low  strain  rates  (<  10~"2/s)  on  universal  testing  machines  or  at  high  strain  rates  (>  102/s)  on 
split  Hopkinson  pressure  bar  (SHPB)  systems.  The  stress-strain  response  of  most  polymeric 
materials  on  the  other  hand  shows  a  pronounced  strain  rate  sensitivity  at  low,  intermediate 
and  high  strain  rates.  At  small  strains,  the  viscoelastic  properties  of  polymers  are  typically 
determined  using  dynamic  mechanical  analysis  (e.g.  McGrumetal.  1997).  The  characterization 
of  the  large  deformation  response  of  polymers  at  low  and  intermediate  strain  rates  of  up  to  10/s 
can  be  performed  on  hydraulic  testing  systems  (e.g.  Yi  et  al.  2006,  Song  et  al.  2007).  As  for 
metals,  conventional  SHPB  systems  are  employed  to  characterize  the  large  deformation  response 
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of  polymeric  materials  at  high  strain  rates.  However,  as  discussed  by  Gray  and  Blumenthal 
(2000),  low  impedance  Hopkinson  bars  are  recommended  when  testing  soft  polymeric  materials 
(e.g.  Zhao  et  al.,  1997,  Chen  et  al.,  1999,  Sharma  et  ah,  2002).  Hoo  Fatt  and  Bekar  (2004) 
developed  a  pulley  system  to  perform  large  strain  tensile  tests  on  rubber  sheets  at  intermediate 
and  high  strain  rates.  Inspired  by  this  work,  Roland  et  al.  (2007)  designed  a  pendulum  impact 
tester  to  study  the  tensile  properties  of  elastomers  at  strain  rates  of  up  to  about  500/s.  In 
both  testing  systems,  the  issues  related  to  the  strain  measurements  under  dynamic  loading 
conditions  are  circumvented  through  the  use  of  digital  image  correlation  (DIC)  based  on  high 
speed  camera  recordings. 

Unlike  for  high  strain  rate  experiments,  the  duration  of  the  experiment  poses  a  major 
challenge  when  using  SHPB  systems  for  intermediate  strain  rate  testing.  The  experiment 
duration  Texp  is  given  by  the  ratio  of  the  strain  emax  at  the  end  of  the  experiment  and  the 
average  strain  rate  e,  Texp  =  £max/e.  In  order  to  avoid  the  superposition  of  waves,  the  maximum 
duration  of  reliable  measurements  is  limited  to  the  input  bar  transit  time.  The  input  bar 
transit  time  is  an  intrinsic  property  of  the  input  bar  and  can  only  be  lengthened  by  increasing 
the  bar  length  or  by  choosing  a  bar  material  of  low  wave  propagation  speed.  In  combination 
with  two  strain  measurements  on  each  Hopkinson  bar,  wave  separation  techniques  may  be 
used  to  overcome  this  limitation  for  elastic  (e.g.  Lundberg  and  Henchoz,  1977,  Yanagihara, 
1978,  Park  and  Zhou,  1999)  and  viscoelastic  bar  systems  (e.g.  Zhao  and  Gary,  1997,  Bacon 
1999,  Casern  et  al.,  2003).  However,  Jacquelin  and  Hamelin  (2001,  2003)  as  well  as  Bussac  et 
al.  (2002)  have  shown  that  so-called  two-point  measurement  wave  separation  techniques  are 
sensitive  to  noise.  This  finding  led  to  the  development  of  a  mathematical  framework  for  an 
advanced  wave  deconvolution  technique  which  is  based  on  redundant  measurements  (Bussac  et 
al.,  2002).  Othman  and  Gary  (2007)  demonstrated  the  applicability  of  this  testing  technique  to 
the  intermediate  strain  rate  testing  of  aluminum  on  a  hydraulic  actuator  driven  SHPB  system. 
Othman  et  al.  (2009)  also  employed  this  technique  when  using  a  0.82m  long  bar  to  measure 
the  axial  forces  in  a  modified  servo-hydraulic  machine.  In  the  present  work,  we  make  use  of  a 
similar  testing  system  as  Othman  and  Gary  (2007)  to  characterize  the  intermediate  strain  rate 
response  of  the  elastomeric  material  polyurea  under  compressive  loading. 

Various  authors  published  experimental  results  on  the  strain  rate  sensitive  response  of  amor- 
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phous  glassy  polymers  (e.g.  Chou  et  al. ,  1973,  Boyce  et  al.,  1988,  Walley  et  al.,  1989,  Cady 
et  al.,  2003,  Siviour  et  al.,  2005,  Mulliken  and  Boyce,  2006,  Mulliken  et  al.,  2006),  crystalline 
glassy  polymers  (e.g.  Chou  et  al.,  1973,  Bordonaro  and  Krempl,  1992,  Cady  et  al.,  2003,  Siviour 
et  al.,  2005,  Khan  and  Farrokh,  2006)  and  elastomers  (e.g.  Gray  et  al.,  1997,  Rao  et  al.,  1997, 
Song  and  Chen,  2003,  2004,  Hoo  Fatt  and  Bekar,  2004,  Shergold  et  al.,  2006,  Roland,  2006)  in¬ 
cluding  polyurea  (e.g.  Amirkhizi  et  al.,  2006,  Roland  et  al.,  2007,  Sarva  et  al.,  2007).  However, 
only  few  experimental  studies  deal  with  the  intermediate  strain  rate  behavior  of  elastomers  at 
large  deformations.  Sarva  et  al.  (2007)  performed  intermediate  strain  rate  compression  tests 
on  polyurea  for  maximum  strains  greater  than  1.0,  but  the  strain  rates  were  only  14  ~  80/s. 
The  same  research  group  also  obtained  test  results  for  a  strain  rate  of  800/s  using  a  very  long 
aluminum  SHPB  system.  Roland  et  al.  (2007)  characterized  the  tensile  behavior  of  polyurea 
over  a  strain  rate  range  of  14  ~  573/s  and  up  to  strains  of  more  than  3.0.  In  the  present  study, 
an  attempt  is  made  to  cover  a  similar  range  of  strain  rates  by  using  the  modified  SHPB  system 
of  Zhao  and  Gary  (1997)  in  combination  with  the  deconvolution  method  of  Bussac  et  al.  (2002) 
to  perform  compression  experiments  on  polyurea. 

This  paper  is  organized  as  follows.  Section  2.2  describes  all  experimental  procedures,  notably 
the  conventional  SHPB  and  the  modified  SHPB  systems.  The  experimental  results  on  polyurea 
are  presented  in  Section  2.3,  followed  by  a  discussion  of  the  limitations  of  the  present  testing 
system  in  Section  2.4. 


2.2  Experimental  Procedures 

Three  different  testing  systems  are  used  to  cover  a  wide  range  of  strain  rates:  a  universal  testing 
machine,  a  conventional  SHPB  system,  and  a  modified  SHPB  system  with  a  hydraulic  actuator. 
Throughout  our  presentation  of  the  experimental  methods,  we  use  the  hat  symbol  to  denote 
the  Fourier  transforms  /  (ui)  =  /  ( t )  e~lujtdt  of  time-dependent  functions  /  ( t ). 

2.2.1  Universal  Testing  Machine 

A  hydraulic  universal  testing  machine  (Model  8800,  Instron,  Canton  MA)  is  used  to  perform 
compression  tests  at  low  and  intermediate  strain  rates  (10-2  ~  10/s).  The  position  of  the 
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vertical  actuator  is  controlled  using  the  software  MAX  (Instron,  Canton).  The  axial  force  F  (t) 
is  measured  using  a  low  profile  load  cell  of  a  maximum  loading  capacity  of  10 kN  (MTS,  Chicago, 
IL)  that  has  been  positioned  at  a  distance  of  25 mm  from  the  specimen.  At  the  same  time,  the 
cross-head  displacement  is  measured  using  an  LVDT  positioned  at  a  distance  of  about  1300mm 
above  the  specimen  (integrated  in  the  actuator  piston).  A  DIC  system  (Vic2D,  Correlated 
Solutions,  Columbia,  SC)  is  employed  to  measure  the  displacements  uul  (t)  and  u0Ut  (t)  of  the 
top  and  bottom  loading  platens,  respectively.  Furthermore,  we  make  use  of  the  DIC  system  to 
quantify  the  Poisson’s  ratio  of  polyurea.  Both  a  thin  polymer  layer  (Teflon)  and  grease  are  used 
to  minimize  the  frictional  forces  at  the  contact  surface  between  the  specimen  and  the  loading 
platens. 


2.2.2  Conventional  SHPB  Systems 

Two  distinct  conventional  SHPB  systems  are  used  in  this  study: 

•  The  first  is  an  aluminum  bar  system  with  a  1203mm  long  striker  bar.  Experiments  of  a 
maximum  duration  of  Texp  =  472 fxs  can  be  performed  on  this  system. 

•  The  second  SHPB  system  is  composed  of  thermoplastic  nylon  bars  with  a  1092mm  long 
striker  bar.  Thus,  a  maximum  duration  of  Texp  =  1255/xs  is  achieved  on  that  system. 

Technical  details  of  these  systems  are  given  in  Table  2.1.  The  gages  for  strain  history 
recordings  are  positioned  near  the  center  of  the  input  bar  and  near  the  output  bar/specimen 
interface  (Table  2.1,  Figure  2-la).  Using  viscoelastic  wave  propagation  theory  (e.g.  Zhao  and 
Gary,  1997),  we  reconstruct  the  incident  and  reflected  waves  based  on  the  input  bar  strain 
gage  recordings  to  estimate  the  force  Fin  (t)  and  displacement  (t)  at  the  input  bar/specimen 
interface.  Analogously,  the  force  Fout  ( t )  and  displacement  uout  (t)  at  the  output  bar/specimen 
interfaces  are  calculated  after  reconstructing  the  transmitted  wave  based  on  the  output  bar 
strain  gage  recording. 

2.2.3  Modified  SHPB  with  Hydraulic  Actuator 

The  total  duration  of  the  loading  pulse  in  an  experiment  on  a  conventional  SHPB  system  is 
limited  by  the  length  of  the  striker  bar.  Thus,  it  is  usually  impossible  to  reach  large  strains  at 
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Figure  2-1:  (a)  Conventional  and  (b)  Modified  SHPB  systems. 


Table  2.1:  Specifications  of  the  conventional  SHPB  sytems 


Aluminum  Bar  System 
Striker  Input  bar  Output  bar 

Striker 

Nylon  Bar  System 

Input  bar  Output  bar 

Length,  L  [m] 

1.203 

2.991 

1.850 

1.092 

3.070 

1.919 

Radius,  R  [m] 

20 

20 

20 

16.5 

20 

20 

Longitudinal  wave 
speed,  co  (cu  =  0) 

[m/s] 

5100 

5100 

5100 

1740 

1740 

1740 

Mass  density,  p  [kg/m3] 

2820 

2820 

2820 

1187 

1162 

1145 

Distance  between  strain 
gauge  and  specimen/bar 
interface,  d  [m] 

- 

1.493 

0.335 

- 

1.537 

0.394 

21 


intermediate  strain  rates  on  conventional  SHPB  systems.  To  overcome  this  key  limitation,  we 
make  use  of  the  modified  SHPB  system  proposed  by  Zhao  and  Gary  (1997).  By  substituting 
the  striker  bar  through  a  hydraulic  actuator,  almost  infinite  loading  pulse  durations  may  be 
achieved.  Having  this  setting,  the  right  end  of  the  output  bar  needs  to  be  fixed  in  space  as 
its  inertia  is  no  longer  sufficient  to  support  the  specimen  (Figure  2- lb).  In  order  to  prevent 
the  failure  of  the  nylon  bars  under  excessive  loads  (it  can  be  difficult  to  stop  the  piston),  a 
fixed  end  support  system  is  designed  such  that  the  bars  are  released  before  elastic  buckling 
occurs.  Note  that  the  wave  superposition  in  the  input  and  output  bars  can  no  longer  be 
avoided  when  the  test  duration  Texp  exceeds  the  transit  time  for  waves  traveling  from  one 
bar  end  to  the  other.  Therefore,  a  wave  separation  technique  is  employed  to  reconstruct  the 
rightward  and  leftward  travelling  waves  in  the  bars  based  on  strain  gage  measurements.  Once 
both  the  rightward  and  leftward  traveling  waves  in  the  bars  are  known,  the  interface  forces  and 
velocities  may  be  calculated  using  the  same  equations  as  those  for  the  input  bar  in  a  conventional 
SHPB  system.  Wave  separation  techniques  in  the  time  domain  are  efficient  for  non-dispersive 
bars  (e.g.  Lundberg  and  Henchoz,  1977),  but  these  require  more  intense  computations  for 
waves  in  dispersive  systems  (e.g.  Bacon,  1999).  Here,  we  adopt  the  frequency  domain  based 
deconvolution  technique  of  Bussac  et  al.  (2002).  In  particular,  displacement  measurements  are 
included  in  addition  to  strain  gage  recordings. 

Suppose  that  a  strain  wave  e(x,t)  in  a  bar  is  composed  of  the  rightward  traveling  wave 
er  ( x,t )  and  the  leftward  traveling  wave  El  ( x,t ).  In  terms  of  Fourier  transforms,  we  have  the 
multiplicative  decomposition  of  the  frequency  and  spatial  dependence, 

e(x,u;)  =£i?(x0,o;)e-^x-IO)+£L(xo,cj)ei«^("-xo)  (2.1) 


where  er  (xo ,  w)  and  el  (xq  ,  w)  denote  the  Fourier  transform  of  the  respective  strain  histories 
at  some  reference  location  xo-  Moreover,  £  (w)  is  the  frequency- dependent  wave  propagation 
coefficient  for  the  respective  bar  system, 


CJ 


£  (w)  =  k  (u>)  +  ia  (w)  =  y  +  ia  (u) 


(2.2) 


with  the  frequency  dependent  wave  number  k(cj),  the  longitudinal  wave  propagation  speed 
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c{iS)  —  uj/k  (to),  and  the  attenuation  coefficient  S  (w). 

The  simplest  approach  to  determine  the  functions  £r  (xq,oj)  and  £/,  ( xq,uj )  is  to  measure  the 
strain  histories  e(x\,t)  and  e(x2,f)  associated  with  the  wave  e(x,t)  at  two  distinct  locations 
x\  and  X2  (on  the  same  bar).  Subsequently,  one  can  solve  the  linear  system  of  equations 


b  =  Ax 


(2.3) 


with  the  unknowns 


the  measurements 


and  the  coefficient  matrix 


x(w) 


£r{x0,u )) 
£l (zq  ,w) 


b(w) 


£  {X\,U>) 
£  (X1,W) 


A(w)  = 


e-if(u)(il-io)  ei£(w)(x  i  —Wo) 
e-i£{u>)(x  2-x0)  ei$(u>)(x  2-x0) 


However,  the  coefficient  matrix  becomes  singular  (det  A  =  0)  if 


(2.4) 


(2.5) 


(2.6) 


£  M  = 


mr 

X2  -  X\ 


(2.7) 


Bussac  et  al.  (2002)  propose  an  integration  method  in  the  complex  domain  to  address  this 
problem.  However,  the  same  authors  have  also  shown  that  the  noise  in  the  recorded  strain  gage 
signals  may  still  lead  to  erroneous  solutions  for  x  (u>).  In  order  to  improve  the  solution  of  Eq. 
(2.3)  under  the  presence  of  measurement  noise,  it  is  useful  to  introduce  redundant  measurements 
including  force,  velocity  or  displacement  measurements.  From  Eq.  (2.1),  the  Fourier  transform 
of  force,  velocity  and  displacement  can  be  expressed  as 


F  (s,  w)  =  E*  (w)  A  [e*  (x0,  u)  e-«M(*-*o) 


u(x, cu)  =  c*  (u)  \-eR(x 0)lj)e-^x-x^ 


+  £l  (x0,u)e^x~^ 
+  £l  (x0,cn)e^(^l 


(2.8) 

(2.9) 
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=  r  (w)  [eR(x0,oj)e-l^x-x^  -  w)  (2.10) 


where 


E*  (u)  =  p 


’  C*(w)  =  ^)  and  l*{u)  =  W) 


(2.11) 


In  the  present  work,  we  perform  only  strain  and  displacement  measurements.  Formally,  we 


e(.xi,w) 


b(w)  = 


e(*Q.w) 

u(xq+i,w) 


(2.12) 


[  u(xQ+R,u)  J 

where  the  subscripts  Q  and  R  represent  the  number  of  measurements  for  strains  and  displace¬ 
ments,  respectively,  at  the  locations  xt  (i  =  1  ,...,Q  +  R)  on  the  bar.  The  corresponding  matrix 
A  (cj)  reads 


o-i£(u)(xi-x0) 


J€(u)(xi-x0) 


a  n 


3-if(oj)(xQ-x0) 


J£(v)(xq-xo) 


l*  (u>)  e-*€H(*Q+l-®o)  _l*  e^(w)(xQ+i-x0) 


(2.13) 


l*  (w)  e-^(w)(x£3+fl-x0)  _£*  e*?M(xQ+K-x o) 

For  redundant  measurements,  the  equation  b  =  Ax  for  the  unknown  x  is  over-determined, 
and  cannot  be  solved  exactly.  Instead,  an  approximate  solution  is  calculated  by  using  the  least 
squares  method  to  minimize  the  scalar  error,  e  =  ||b  —  Ax||2  =  (b  —  Ax)^  (b  —  Ax).  Thus, 
the  approximate  solution  x  minimizing  the  error  must  satisfy  the  equation  AHb  =  AflAx.  As 
long  as  the  columns  of  A  are  linearly  independent,  the  matrix  AH  A  is  positive  definite  (e.g. 
Strang,  1985)  and  the  unknown  x  can  be  determined  as 


x  =  [AH  A]  1  AHb 


(2.14) 
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where  the  Hermitian  AH  (complex  conjugate  and  transpose  of  A)  corresponds  to  the  transpose 
of  A  if  A  is  real  (e.g.  Magnus  and  Neudecker,  1988).  Note  that  a  least  squares  solution  of 
similar  form  has  been  presented  by  Hillstrom  et  al.  (2000)  in  the  context  of  complex  modulus 
identification  based  on  redundant  strain  measurements. 

In  order  to  rule  out  the  linear  dependence  of  the  columns  of  A,  we  modify  the  propagation 
coefficient  £  (w)  artificially.  In  other  words,  when  calculating  A,  £  (w)  is  substituted  by  the 
modified  propagation  coefficient  £  (to) 

(2.15) 

where  77  is  a  very  small,  but  otherwise  arbitrary,  negative  number;  throughout  our  analysis, 
we  used  fj  =  — 10-7.  The  modified  propagation  coefficient  £  (co)  corresponds  to  the  linear 
perturbation  of  the  propagation  coefficient  £  (to)  in  the  complex  frequency  domain.  As  a  result, 
the  propagation  coefficient  is  always  complex- valued,  and  thus  the  singularity  condition  in  Eq. 
(2.7)  can  no  longer  be  satisfied.  Note  that  fj  introduces  a  very  small  artificial  attenuation 

~.  ,  f  .  .  _da(u>)]  .  T  ,  .  _d,Kluj)  1 

£  (u>)  =  — - -  +  1  a(io)+r]— -  (2-16) 

CLUJ  |_  CLUJ 

As  a  result,  even  purely  elastic  materials  ( a  (w)  =  0)  exhibit  some  artificial  attenuation  (i.e. 
non-zero  imaginary  value)  which  ensures  the  causality  of  the  waves  propagating  in  a  bar  (e.g. 
Bacon,  1999). 

In  the  present  study  we  make  use  of  a  modified  SHPB  system  with  nylon  input  and  output 
bars.  Table  2.2  summarizes  the  technical  specifications  of  the  testing  system.  Each  bar  is 
equipped  with  three  strain  gages  and  a  high  contrast  grid  for  optical  displacement  measurements 
(Model  100H,  Zimmer,  Germany).  After  using  Eq.  (2.14)  to  determine  the  leftward  and  the 
rightward  traveling  waves  e™  (xl0 (l,t)  and  {x o\i),  the  displacement  Uin(t)  and  the  force 
Fin  ( t )  at  the  input  bar/specimen  interface  are: 

Uin  M  -  r  M  [eft  (xjn,  w)  e^o"  -  ef  e-*(“Kn]  (2.17) 

Fm  (oj)  =  E*  (w)  A  [eft  (xft,  u)  +  eft  (xft,  u)  (2.18) 
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Table  2.2:  Specifications  of  the  modified  nylon  SHPB  sytem 


Input  bar 

Output  bar 

Length,  L  [m] 

3.123 

3.045 

Radius,  R  [m] 

20 

20 

Longitudinal  wave 
speed,  Co  (w  =  0) 

[■ m/s ] 

1740 

1740 

Mass  density,  p  [kg/m3] 

1150 

1150 

0.61 

0.825 

Distance  between  strain 
gauges  and  specimen/bar 
interface,  dsg  [m] 

1.515 

1.523 

2.623 

2.582 

Distance  between  displacement 
sensor  and  specimen/bar 
interface,  ddm  [rri] 

0.953 

2.183 

Analogously,  the  displacement  uout  (uj)  and  the  corresponding  force  Fout  (to)  at  the  output 
bar/specimen  interface  are  determined  from  the  output  bar  measurements. 

2.2.4  Determination  of  the  Stress-Strain  Curves 

The  time  histories  of  the  displacements  and  forces  at  the  specimen  boundaries,  Uin  ( t ),  uout  ( t ), 

Fin  ( t )  and  Fout  (f)  are  obtained  from  applying  the  inverse  Fourier  transform  /  (t)  =  ^  /  (w)  elUJtdcu 

to  Uin  (w),  u0ut  (w),  Fin  (w)  and  Fout  (w),  respectively.  The  input  force  Fjn  (t)  is  considered  as 
a  redundant  measurement;  it  is  used  to  verify  the  condition  of  quasi-static  equilibrium  of  a 
dynamically  loaded  specimen 

Fin  ( t )  -  Font  ( t )  ^  0  (2.19) 
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The  spatial  average  of  the  logarithmic  axial  strain  within  the  specimen  reads 


£  (t)  =  In  (l  +  (2.20) 

where  Lq  denotes  the  initial  length  of  the  specimen.  Using  the  output  force  measurement,  we 
calculate  the  true  stress 

a  (t)  =  I'"u[  exp  [2 us  (f)]  (2-21) 

Aq 

where  Aq  is  the  initial  cross-sectional  area,  and  v  is  the  elastic  Poisson’s  ratio.  In  the  present 
work,  it  is  assumed  that  the  Poisson’s  ratio  is  constant,  i.e.  it  depends  neither  on  the  strain 
nor  on  the  strain  rate.  The  final  true  stress-strain  curve  is  then  found  from  the  combination  of 
the  stress  and  strain  history  functions 


a  (e)  —  a  (t)  o  e  (t)  (2.22) 

2.3  Experimental  Results 

All  experiments  are  performed  on  polyurea.  This  rubber-like  material  has  a  mass  density  of 
1.0 gjcrv?  and  an  elastic  modulus  of  about  lOOMPa. 

2.3.1  Experiments  using  the  Universal  Testing  Machine 

Representative  stress-strain  curves  for  true  strain  rates  of  up  to  10/s  are  determined  from 
experiments  on  the  universal  testing  machine  using  cylindrical  specimens  of  diameter  Do  — 
10mm  and  length  Lq  =  10mm.  All  experiments  are  carried  out  under  displacement  control  up 
to  a  maximum  true  compressive  strain  of  1.0  (which  corresponds  to  an  engineering  compressive 
strain  of  0.63).  In  order  to  achieve  a  constant  true  strain  rate  £o,  a  velocity  profile  Uin  (*)  of 
exponential  shape  is  applied  to  the  top  of  the  specimen 

Uin  {t)  =  ~Lq£q  exp  (e0t)  (2.23) 

The  Poisson’s  ratio  is  determined  from  the  experiments  performed  at  true  compressive  strain 
rates  of  up  to  1/s.  Based  on  the  DIC  measurements  of  the  specimen  diameter  D  =  D(t),  we 
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Logarithmic  Axial  Strain  [-] 


Figure  2-2:  Identification  of  Poisson’s  ratio  from  the  linear  relationship  between  the  logarithmic 
radial  and  axial  strains 

calculate  the  logarithmic  radial  strain  eT , 


where  Dq  denotes  the  initial  specimen  diameter.  The  experimental  data  depicted  in  Figure  2-2 
shows  the  linear  relationship  between  the  logarithmic  radial  strain  and  the  logarithmic  axial 
strain.  Upon  evaluation  of  the  slope,  we  find  a  Poisson’s  ratio  of  u  —  0.448. 

The  data  acquisition  rate  of  the  DIC  system  is  limited  to  about  1Hz.  Thus,  we  only  use  the 
DIC  system  for  the  slowest  experiments  and  make  use  of  the  actuator  position  measurement 
(LVDT)  to  determine  the  effective  axial  displacement  at  higher  strain  rates.  The  comparison 
of  the  LVDT  readout  with  the  DIC  measurement  yields  an  overall  stiffness  of  the  testing  frame 
of  about  100 kN/mm.  The  measured  true  stress-strain  curves  are  shown  in  Figure  2-3a  for 
true  strain  rates  of  about  10_2/s,  10-1/s>  10°/s,  and  lO1/^  The  corresponding  true  strain  rate 
versus  true  strain  curves  are  depicted  in  Figure  2-3b.  For  the  slowest  experiment  (e o  —  10~2/s), 
the  slope  of  the  stress-strain  curve  (Figure  2-3a)  decreases  significantly  at  a  stress  of  about  0.1; 
subsequently,  the  stress-strain  curve  changes  its  shape  from  concave  to  convex  at  an  axial 
strain  of  about  0.3.  Due  to  the  characteristic  rubber  chain  locking  behavior,  the  stress  level 
increases  monotonically  throughout  the  entire  experiment  from  6 MPa  at  e  =  0.1  to  13.5 MPa 
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Figure  2-3:  Test  results  from  the  universal  testing  machine:  (a)  True  stress-strain  curves,  (b) 
True  strain  versus  true  strain  curves 

at  e  =  1.0.  For  the  next  higher  strain  rate  (e o  —  10_1/s),  the  overall  stress  level  is  about  12% 
higher.  Similarly,  the  shape  of  the  stress-strain  curve  is  preserved  for  strain  rates  of  10°/s  and 
lO1/^,  but  the  stress  level  increases  by  35%  and  65%,  respectively,  compared  to  that  at  10_2/s. 


2.3.2  Experiments  Using  Conventional  SHPB  Systems 

Appendix  B  outlines  the  identification  of  the  frequency  dependant  coefficients  c  (tu)  and  a  (oj) 
for  both  the  aluminum  and  nylon  bars.  The  results  axe  presented  in  Figure  2-4  together  with 
the  Pochhammer-Chree  solution  (e.g.  Graff,  1975).  These  experimentally  obtained  coefficients 
are  used  throughout  our  analysis  of  the  waves  in  both  the  conventional  and  the  modified  SHPB 
systems. 

Aluminum  Bar  System 

Experiments  at  high  strain  rates  are  performed  on  the  conventional  aluminum  SHPB  system. 
Cylindrical  polyurea  specimens  with  Do  =  20 mm  and  Lq  =  5 mm  are  used  on  the  aluminum 
system.  Average  strain  rates  of  e  ~  3700/s  and  i  —  2300/s  are  achieved  at  striker  velocities 
of  13 m/s  and  9 m/s,  respectively.  To  verify  the  quasi-static  equilibrium  throughout  the  ex¬ 
periments,  both  the  input  and  output  force  are  depicted  in  Figures  2-5a  and  2-5b.  The  poor 
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Figure  2-4:  Propagation  coefficient  of  aluminum  bar  (first  row)  and  nylon  bar  (second  row). 
(a)-(c)  Longitudinal  wave  speed  and  (b)-(d)  Attenuation  coefficient. 
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agreement  of  the  force  measurements  for  3700/s  may  be  read  as  lack  of  equilibrium  (e.g.  Aloui 
et  al.,  2008).  However,  for  the  present  experiments,  this  observation  is  attributed  to  the  low 
signal-to-noise  ratio  for  the  input  force  measurements.  Due  to  the  pronounced  mismatch  be¬ 
tween  the  force  amplitude  of  the  incident  wave  (e.g.  Finc  —  AalP alcal^ 'str /2  —  120kN  for 
13 m/s)  and  the  specimen  resistance  (Fln  —  Aspcaspc  ~  15kN  at  e  =  0.5),  most  of  the  inci¬ 
dent  wave  is  reflected  at  the  input  bar /specimen  interface  which  ultimately  results  in  a  poor 
input  force  measurement  (see  e.g.  Grolleau  et  al.  (2008)  for  details  on  the  force  measurement 
accuracy).  The  incident  wave  exhibits  some  Pochhammer-Chree  oscillations  due  to  the  lateral 
inertia  of  the  40mm  diameter  aluminum  bars.  Consequently,  we  observe  some  non-monotonic 
behavior  in  the  stress-strain  curves  for  3700/s  and  2300/s  in  Figure  2-6a.  The  overall  shape  of 
the  curve  is  very  similar  to  that  for  static  loading,  but  the  stress  level  is  almost  three  to  four 
times  higher. 

Nylon  Bar  System 

Another  set  of  high  strain  rate  experiments  (1200/s  and  700/s)  are  performed  using  smaller 
diameter  specimens  (Do  =  10 mm,  Lq  =  5 mm  )  on  the  nylon  bar  SHPB  system.  Recall  that  the 
main  reason  for  changing  from  aluminum  to  nylon  bars  is  to  increase  the  maximum  duration  of 
the  experiments  from  Texp  =  472 fis  to  1255/rs.  At  the  same  time,  the  use  of  nylon  significantly 
reduces  the  impedance  mismatch  between  the  bars  and  the  polyurea  specimen.  This  improves 
the  force  measurement  accuracy,  notably,  that  of  the  input  force.  Striker  velocities  of  8 m/s  and 
6 m/s  are  needed  to  obtain  an  average  strain  rate  of  e  cs  1200/s  and  t  —  700/s,  respectively. 
Higher  striker  velocities  would  cause  inelastic  deformation  in  the  bars  upon  striker  impact.  On 
the  other  hand,  for  a  maximum  loading  duration  of  1255/us,  lower  striker  velocities  would  not 
achieve  the  desired  maximum  true  compressive  strain  of  e  =  1.0. 

There  are  less  signal  oscillations  in  the  nylon  than  in  the  aluminum  system  because  of 
its  higher  signal-to-noise  ratio.  Furthermore,  due  to  the  lower  striker  velocities  and  the  wave 
attenuation  in  the  nylon  input  bar,  there  are  less  severe  Pochhammer-Chree  oscillations  in 
the  incident  wave  signal  as  compared  to  the  aluminum  system  (see  Figure  2-55).  Therefore, 
relatively  smooth  stress-strain  curves  are  obtained  from  the  dynamic  experiments  on  the  nylon 
bar  system  (Figure  2-6a).  The  exact  evolution  of  the  true  strain  rates  as  a  function  of  the 
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Figure  2-5:  Comparison  of  forces  between  input  bar  and  output  bar  from  SHPB  tests  from 
aluminum  bar  tests:  (a)  3700/s,  (b)  2300/s;  and  nylon  bar  tests:  (c)  1200/s,  (d)  700/s 
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Figure  2-6:  (a)  True  stress-strain  curves,  (b)  True  strain  rate  versus  true  strain  curves. 

true  strain  are  shown  in  Figure  2-6b.  Unlike  for  the  experiments  on  the  aluminum  bar  system, 
the  true  strain  rate  is  no  longer  increasing  in  a  monotonic  manner.  This  is  due  to  the  lower 
amplitude  of  the  incident  force  (e.g.  for  F{nc  —  AalP ALcALYstr /2  —  7.6 kN  for  e  ~  700/s)  which 
is  now  of  the  same  order  of  magnitude  as  the  specimen  resistance  (Fin  —  Aspctjspc  —  2.5 kN). 
Thus,  as  the  specimen  resistance  increases  throughout  the  experiment,  the  magnitude  of  the 
reflected  wave  decreases;  as  a  result,  despite  the  logarithmic  strain  definition,  the  engineering 
strain  rate  no  longer  increases  due  to  the  decreasing  interface  velocity. 

2.3.3  Experiments  Using  the  Modified  SHPB  Systems 

Experiments  are  performed  on  the  modified  nylon  SHPB  system  using  the  hydraulic  actuator 
in  an  open  mode,  which  is  different  from  the  conventional  closed  loop  mode  of  servo  hydraulic 
testing  machines.  In  this  open  loop  mode,  the  user  can  preset  the  position  of  the  inlet  servo 
valve.  Furthermore,  the  initial  pressure  of  the  in-flowing  fluid  may  be  controlled.  However, 
the  user  has  no  active  control  of  the  actuator  velocity  throughout  the  experiment.  Actuator 
piston  velocities  of  up  to  5 m/s  may  be  achieved  in  this  mode  of  operation.  Here,  we  perform 
experiments  at  Am/ s,  1  m/s,  0.5 m/s,  and  0.1  m/s  which  resulted  in  average  compressive  strain 
rates  of  about  1000/s,  110/s,  36/s  and  10/s. 

Three  strain  gages  and  one  displacement  measurement  are  taken  into  account  (per  bar)  to 
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reconstruct  the  waves  in  either  bar  using  the  above  deconvolution  technique1 .  The  comparison  of 
the  measured  input  and  output  force  histories  confirms  the  quasi-static  equilibrium  for  110/s, 
36/s  and  10/s  (Figure  2-7).  The  differences  between  the  input  and  output  force  for  1000/s 
are  associated  with  the  poor  quality  of  the  deconvolution  based  estimate  of  the  input  force; 
the  accuracy  of  the  optical  displacement  measurement  system  decreases  substantially  at  high 
loading  velocities  leading  to  severe  oscillations  in  the  input  force  history.  However,  considering 
that  the  higher  velocity  cases  (Figure  2-5)  show  the  good  force  agreement,  the  quasi-static 
equilibrium  can  also  be  assumed  for  the  strain  rate  of  1000/s.  A  significant  force  drop  is  found 
at  f  ~  5 ms,  20 ms,  and  60ms  for  the  strain  rates  of  110/s,  36/s  and  10/s,  respectively.  This 
force  drop  is  due  to  the  premature  partial  failure  of  the  fixed  end  support  of  the  output  bar 
that  causes  a  short  unloading-reloading  cycle.  The  same  force  drops  are  also  found  in  the  stress 
versus  strain  curve  (Figure  2-8a)  at  strains  of  0.55,  0.60  and  0.25  for  strain  rates  of  110/s,  36/s 
and  10/s,  respectively. 

2.3.4  Comment  on  the  Signal  Oscillations 

There  are  two  characteristic  time  scales  associated  with  the  experiments  on  the  modified  SHPB 
system.  The  short  time  scale  corresponds  to  the  round  trip  time  of  an  elastic  wave  traveling 
through  the  specimen,  Tspc  cs  11/rs.  The  second  time  scale  is  much  longer;  it  is  associated  with 
the  round  trip  of  an  elastic  wave  traveling  through  the  input  bar,  T*n  ~  3600^ts.  The  experiment 
at  an  average  strain  rate  of  1000/s  remains  unaffected  by  the  large  time  scale  as  the  total 
duration  of  the  experiment  ( Texp  cs  10 3/xs)  is  still  shorter  than  T)n  ~  3600 /is.  However,  already 
at  a  strain  rate  of  110/s,  the  duration  of  the  experiment  ( Texp  ~  104/xs)  exceeds  Tln  ~  3600/r.s. 
As  a  result,  the  shape  and  amplitude  of  the  incident  wave  is  not  only  determined  by  the  velocity 
of  the  hydraulic  actuator,  but  also  by  the  leftward  travelling  wave  that  has  been  reflected  by 
the  specimen.  Consequently,  the  incident  wave  changes  with  a  periodicity  of  T\n.  In  the  present 
experiments,  the  first  reflected  wave  is  a  tensile  wave  which  reduces  the  initial  magnitude  of 
the  compressive  incident  wave.  Hence,  the  rate  of  loading  decreases  before  the  rate  of  loading 
increases  again  after  the  next  period  of  Tm.  Therefore,  this  abrupt  change  of  the  loading  rate 

'The  only  exception  is  the  input  bar  in  the  experiment  at  110/s  where  one  of  the  three  strain  gauge  signals 
was  not  properly  recorded.  Therefore,  only  two  strain  gages  and  one  displacement  measurement  were  taken  into 
account  for  that  experiment. 
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Figure  2-7:  Comparison  of  forces  between  input  bar  and  output  bar  from  the  modified  SHPB 
tests:  (a)  1000/s,  (b)  110/s,  (c)  36/s  and  (d)  10/s 
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Figure  2-8:  Test  results  from  the  modified  SHPB  system:  (a)  True  stress-strain  curves,  (b)  True 
strain  rate  versus  true  strain  curves. 

has  a  periodicity  of  T)n .  The  corresponding  strain  rate  versus  strain  curve  shows  a  pronounced 
decrease  in  strain  rate;  since  the  strain  increases  only  little  during  a  period  of  reduced  loading 
rate,  we  observe  sharp  drops  in  the  strain  rate  versus  total  strain  curve.  For  lower  average 
strain  rates,  this  number  of  strain  rate  drops  increases  further.  Formally,  we  may  write 

Save 


where  n  is  the  number  of  the  expected  drops  in  strain  rate  associated  with  the  wave  reflections 
in  the  input  bar.  This  number  is  2,  7,  and  27  for  the  experiments  at  average  strain  rates  of 
110/s,  36/s  and  10/s.  In  the  limiting  case  of  static  loading,  we  have  n  — ►  oo  which  ultimately 
results  in  a  constant  strain-rate  versus  strain  curve.  In  addition  to  loading  velocity  changes 
associated  with  wave  reflections  in  the  input  bar,  our  experimental  results  are  also  affected  by 
other  sources  of  vibrations.  These  include  wave  reflections  within  the  output  bar  as  well  as 
vibrations  in  the  fixed  end  support  system  and  the  hydraulic  actuator.  Therefore,  the  exact 
identification  of  all  strain  rate  drops  in  Figure  2-8b  according  to  Eq.  (2.25)  has  been  omitted. 
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Figure  2-9:  Comparison  of  the  results  from  the  modified  SHPB  with  those  from  other  two 
testing  methods:  (a)  True  stress-strain  curves,  (b)  True  strain  rate  versus  true  strain  curves. 

2.4  Discussion 

2.4.1  Experimental  Results 

To  validate  our  experimental  data,  we  first  checked  the  consistency  among  the  results  obtained 
from  different  testing  methods.  Figure  2-9  shows  selected  stress-strain  curves  obtained  from  the 
modified  SHPB  system  (dashed  lines)  next  to  the  results  from  the  conventional  SHPB  (red  solid 
line)  and  the  universal  testing  machine  (blue  solid  line).  For  1000/s,  the  modified  SHPB  result 
shows  reasonably  good  agreement  with  the  conventional  SHPB  curve  for  1200/s.  Analogously, 
for  the  average  strain  rate  of  10/s,  the  stress-strain  curve  obtained  from  the  modified  SHPB 
test  corresponds  well  to  that  obtained  from  the  test  on  the  universal  testing  machine.  Recall 
that  the  perturbation  of  the  stress-strain  curve  for  the  modified  SHPB  system  at  about  0.25 
strain  is  due  to  the  partial  premature  failure  of  the  output  bar  end  support  system.  The  stress 
level  from  the  modified  SHPB  is  slightly  higher  after  partial  support  failure  which  is  attributed 
to  differences  in  the  strain  rate. 

The  data  in  Figure  2- 10a  show  the  stress  as  a  function  of  the  strain  rate  for  different  levels 
of  strain:  0.1,  0.5  and  0.9.  The  effect  of  strain  rate  is  more  pronounced  at  large  strains.  For 
instance,  at  a  strain  of  0.1,  the  stress  level  increases  by  317%  when  increasing  the  strain  rate 
from  e  ~  10~2/s  to  £  ~  3700/s  (increase  6 kN  from  to  19/cAQ;  at  a  strain  of  0.9,  the  stress  level 
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Figure  2-10:  True  stress  as  a  function  of  the  strain  rate  at  selected  strain  levels:  (a)  Results  of 
the  current  study  and  a  fit  of  Eq.  (2.26)  to  the  results  from  the  present  study,  (b)  Comparison 
of  the  results  with  previous  studies. 


increases  by  360%  over  the  same  range  of  strain  rates  (from  12.5 kN  at  k  ~  10_2/s  to  AbkN  at 
k  ~  3700/s).  Using  the  separable  form  of  strain  and  strain-rate  effect,  in  the  same  figure,  the 
following  nonlinear  empirical  function  has  been  fitted  to  the  experimental  data: 


cr(e,£) 


/(e)  5(e) 


/(e) 


1  + A 


(2.26) 


where  /  (e)  is  the  asymptotic  stress-strain  curve  at  infinitely  slow  loading  conditions  and  g  (e) 
represents  the  contribution  of  strain  rate  sensitivity.  Here,  the  reference  strain  rate  of  kref  = 
10~2/s  is  chosen  to  be  the  lowest  strain  rate  in  the  present  experiments.  In  addition,  A  =  0.0899 
is  responsible  for  the  amplitude  of  the  strain  rate  contribution,  and  B  —  0.249  represents  the 
sensitivity  of  the  strain  rate. 


2.4.2  Intermediate  Strain  Rate  Testing  Systems 

The  present  experimental  study  confirms  the  high  strain  rate  sensitivity  of  the  polyurea  material 
which  has  been  reported  in  earlier  studies.  Roland  et  al.  (2007)  performed  a  series  of  tensile 
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tests  using  a  custom-made  pulley  system  in  a  drop  tower  to  perform  uniaxial  tensile  tests  at 
intermediate  and  high  strain  rates.  Sarva  et  al.  (2007)  used  an  enhanced  universal  testing 
machine  to  perform  compression  experiments  at  strain  rates  of  up  to  80/s  while  an  aluminum 
SHPB  system  with  a  striker  bar  length  of  3 m  has  been  used  to  perform  experiments  at  strain 
rates  above  800/s.  The  comparison  of  the  present  experimental  data  with  the  results  of  Sarva 
et  al.  (2007)  and  Roland  et  al.  (2007)  confirms  the  validity  of  the  measurements  with  the 
modified  SHPB  system  (Figure  2-10b). 

The  implementation  of  the  deconvolution  technique  by  Bussac  et  al.  (2002)  leads  to  a  stable 
algorithm  that  is  convenient  to  use  for  the  reconstruction  of  dispersive  waves  in  bars  based 
on  redundant  measurements.  Thus,  the  theoretical  limitation  of  the  duration  of  experiments 
on  SHPB  systems  is  successfully  overcome.  In  combination  with  a  hydraulic  actuator,  the 
entire  range  of  low  to  high  strain  rates  could  be  covered  using  a  single  testing  system.  The 
comparison  with  conventional  SHPB  experiments  at  high  strain  rates  and  universal  testing 
machine  experiments  at  low  strain  rates  has  confirmed  the  validity  of  the  modified  SHPB 
technique.  However,  there  are  still  two  difficulties  associated  with  our  modified  SHPB  system 
which  need  to  be  addressed  in  the  future: 

•  Displacement  and/or  velocity  measurement  accuracy:  the  accuracy  of  the  deconvolution 
technique  relies  heavily  on  accurate  displacement  measurements  (in  particular  at  low 
strain  rates).  The  present  optical  technique  provided  good  results  for  loading  velocities 
of  up  to  0.5?n/s,  but  significant  errors  became  visible  at  larger  loading  velocities. 

•  Quality  of  the  loading  pulse  at  specimen  interfaces.  In  order  to  achieve  approximately  con¬ 
stant  strain  rates,  the  ideal  loading  pulse  should  be  such  that  the  bar/specimen  interfaces 
move  at  constant  velocities. 

The  first  difficulty  may  be  resolved  through  the  use  of  improved  measurement  equipment. 
Alternatively,  the  deconvolution  technique  for  high  loading  velocities  may  also  be  applied  using 
strain  gage  measurements  only.  However,  it  is  very  challenging  to  overcome  the  second  difficulty. 
As  an  intermediate  strain  rate  experiment  takes  much  longer  than  a  wave  round  trip  in  the 
input  bar,  the  input  bar/specimen  interface  velocity  is  not  constant  even  if  the  hydraulic  piston 
moves  at  a  constant  velocity.  Simple  wave  analysis  shows  that  a  period  of  high  velocity  loading 
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is  followed  by  a  period  of  loading  at  a  lower  rate;  the  length  of  each  period  corresponds  to  the 
round  trip  time  for  a  wave  travelling  in  the  input  bar.  The  same  holds  true  for  the  output 
bar/specimen  interface  velocity  which  is  affected  by  the  round  trip  time  in  the  output  bar. 
Consequently,  the  strain  rate  in  our  intermediate  strain  rate  experiments  was  not  constant. 

For  the  desired  maximum  true  compressive  strain  of  e  =  1.0,  the  total  duration  of  an 
intermediate  strain  experiment  at  50/s  is  Texp  —  20 ms  -  irrespective  of  the  specimen  geometry. 
Conceptually,  there  exist  several  solutions  to  this  problem: 

•  Conventional  nylon  SHPB  system  with  a  striker  bar  length  of  1740  x  0.02  x  0.5  =  17.4m 
along  with  a  35m  long  input  bar  and  a  17.5m  long  output  bar.  In  this  configuration,  all 
strain  gages  can  be  positioned  such  that  the  rightward  and  leftward  traveling  waves  do 
not  superpose  at  the  strain  gage  locations. 

•  Conventional  nylon  SHPB  system  with  a  17.5?n  long  striker  bar  and  17.5m  long  input 
and  output  bars.  In  this  case,  a  deconvolution  method  needs  to  be  used  to  reconstruct 
the  waves  in  the  input  and  output  bars.  However,  the  input  bar  is  still  sufficiently  long 
to  guarantee  that  the  round  trip  time  is  greater  than  the  duration  of  the  experiment. 

•  Hydraulic  nylon  SHPB  system  with  17.5m  long  input  and  output  bars.  Based  on  the 
assumption  that  the  hydraulic  piston  moves  at  a  constant  velocity,  this  system  will  provide 
the  same  capabilities  as  the  previous  system. 

As  an  alternative  to  very  long  input  and  output  bars,  one  may  chose  the  opposite  strategy. 
Note  that  the  magnitude  of  the  oscillations  is  proportional  to  the  change  in  force  level  in 
the  specimen  over  the  time  Tin.  Thus,  the  shorter  the  input  bar,  the  smaller  the  oscillation 
magnitude.  One  could  therefore  envision  very  short  (e.g.  <  0.5m)  small  diameter  input  and 
output  bars.  In  this  case,  we  have  Tj„  =  0.57 ms  and  hence  Tin  <C  Texp. 

However,  since  the  modified  SHPB  system  requires  two  displacement  measurement  sensors 
(notably  for  low  strain  rate  experiments),  one  can  also  use  these  sensors  to  measure  directly 
the  displacements  of  the  respective  bar/specimen  interfaces.  Hence  the  strain  history  can  be 
measured  without  using  the  deconvolution  algorithm.  The  bars  would  therefore  only  serve  as 
load  cell  to  measure  the  force  history.  Unless  the  quasi-static  equilibrium  needs  to  be  verified 
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experimentally,  a  single  force  measurement  is  satisfactory.  Moreover,  it  may  be  worth  consider¬ 
ing  a  piezoelectric  sensor  to  measure  the  force,  thereby  completely  eliminating  the  use  of  bars 
to  perform  the  experiments  at  low,  intermediate  and  high  strain  rates.  As  our  hydraulic  piston 
cannot  provide  a  constant  loading  velocity  above  0.5 m/s,  a  striker  bar  may  also  be  used  to  load 
the  specimen.  The  only  unknown  which  is  left  in  this  system  is  the  realization  of  the  “fixed” 
boundary  condition.  Further  research  needs  to  be  carried  out  to  design  a  support  point  that 
does  not  introduce  spurious  oscillation  into  the  testing  system. 

2.5  Conclusion 

The  modified  SHPB  system  of  Zhao  and  Gary  (1997)  has  been  used  to  perform  compression  tests 
on  polyurea  at  low,  intermediate  and  high  strain  rates.  It  is  composed  of  nylon  input  and  output 
bars,  while  the  striker  bar  is  substituted  by  a  hydraulic  actuator.  Using  the  deconvolution 
technique  by  Bussac  et  al.  (2002),  the  time  limitation  of  conventional  SHPB  systems  may  be 
overcome,  thereby  enabling  the  use  of  the  modified  SHPB  system  for  low  and  intermediate  strain 
rate  experiments  of  long  duration.  The  experiments  confirm  the  known  strain  rate  sensitivity 
of  polyurea.  The  measured  stress  levels  correspond  well  to  earlier  results  which  have  been 
obtained  from  tests  on  conventional  SHPB  systems  with  very  long  bars.  Although  the  intrinsic 
time  limitation  of  SHPB  systems  could  be  overcome,  this  study  also  shows  that  it  is  still  not 
possible  to  perform  experiments  at  reasonably  constant  strain  rates  with  this  technique.  This  is 
due  to  the  finite  length  of  the  input  and  output  bars  which  causes  a  periodic  change  in  loading 
velocity.  It  is  shown  that  intermediate  strain  rate  SHPB  experiments  require  either  very  long 
bars  (>  20 m)  or  very  short  bars  (<  0.5m)  in  order  to  achieve  an  approximately  constant  strain 
rate  throughout  the  entire  experiment. 
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Chapter  3 


Constitutive  Modeling 


3.1  Introduction 

Polyurea  is  used  to  mitigate  structural  damage  during  impact  loading  because  of  its  good 
damping  performance.  In  addition,  it  is  utilized  by  various  industries  because  of  its  fast  setting 
time  as  well  as  its  good  chemical  and  fire  resistance.  Polyurea  has  found  applications  in  army 
vehicles  for  blast  protection  because  of  its  high  toughness-to-density  ratio  at  high  strain  rates. 
Polyurea  shows  a  highly  nonlinear  viscoelastic  behavior  at  finite  strains  (e.g.  Amirkhizi  et 
al.,  2006,  Bogoslovov  and  Roland,  2007,  Roland  et  al.,  2007,  Shim  and  Mohr,  2009a).  The 
mechanical  properties  of  linearly  viscoelastic  materials  may  be  described  by  the  relaxation 
modulus  (or  creep  compliance)  which  is  independent  of  strain  magnitude.  However,  nonlinear 
viscoelasticity  is  characterized  by  a  decrease  (or  increase)  of  the  relaxation  modulus  (or  creep 
compliance)  with  increasing  strain  or  decreasing  stress  (e.g.  Brinson  and  Brinson,  2008). 

Most  finite  viscoelasticity  models  of  elastomers  are  formulated  using  either  (1)  the  so-called 
hereditary  integral  approach  or  (2)  the  framework  of  multiplicative  decomposition  of  the  de¬ 
formation  gradient.  Motivated  by  linear  viscoelastic  models,  hereditary  integral  models  are 
formulated  in  terms  of  relaxation  or  memory  functions  (e.g.  Lockett,  1972).  Widely  used  single 
integral  theories  are  the  theory  of  Finite  Linear  Viscoelasticity  (Coleman  and  Noll,  1961)  and 
the  BKZ  theory  (Bernstein  et  al.,  1963);  both  make  use  of  several  relaxation/memory  functions. 
Significant  efforts  have  been  made  to  improve  these  theories  and  to  reduce  the  number  of  re¬ 
quired  material  parameters  (e.g.  Lianis,  1963,  McGuirt  and  Lianis,  1970,  Leonov,  1976,  Johnson 
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et  al.,  1994,  Haupt  and  Lion,  2002).  Nonlinear  viscoelastic  behavior  is  often  considered  as  the 
superposition  of  a  rate-independent  nonlinear  elastic  response  (so-called  equilibrium  part)  and  a 
viscosity-induced  overstress  contribution  which  is  described  through  fading  memory  functions. 
Most  nonlinear  viscoelastic  constitutive  models  have  been  experimentally  validated  at  very  low 
strain-rates  of  0.1/s  or  less.  Only  few  papers  deal  with  the  nonlinear  viscoelastic  behavior  of 
elastomers  at  intermediate  and  high  strain-rates  (1/s  to  1000/s).  Using  the  hereditary  integral 
approach,  Yang  et  al.  (2000)  and  Shim  et  al.  (2004)  proposed  a  phenomenological  constitutive 
model  to  predict  the  behavior  of  silicon  rubber  at  high  strain  rates  (900/s  to  3000/s).  Hoo  Fatt 
and  Ouyang  (2007)  adopted  the  integral  approach  to  model  the  response  of  butadiene  rubber  at 
strain  rate  ranging  from  76/s  to  450/s.  The  validity  of  most  hereditary  integral  approach  based 
models  is  limited  to  a  narrow  range  of  strain  rates  due  to  the  use  of  only  one  relaxation  time 
period  (Yang  et  al.,  2000  and  Shim  et  al.,  2004)  or  a  constant  memory  function  (Hoo  Fatt  and 
Ouyang,  2007).  In  general,  the  hereditary  integral  approach  is  very  useful  in  describing  finite 
viscoelastic  behavior,  but  its  successful  application  to  the  real  test  data  depends  strongly  on 
the  effectiveness  of  the  rather  complex  relaxation  or  memory  function  calibration  procedures. 

Compared  to  hereditary  integral  models,  the  multiplicative  decomposition  of  the  deforma¬ 
tion  gradient  typically  leads  to  models  with  material  parameters  that  can  be  easily  identified 
from  experiments.  The  concept  of  the  multiplicative  decomposition  of  the  deformation  gra¬ 
dient  (Kroner,  1960,  and  Lee,  1969)  was  initially  applied  to  finite  viscoelasticity  by  Sidoroff 
(1974)  and  further  explored  by  others  (e.g.  Lubliner,  1985,  Le  Tallec,  1993,  Reese  and  Govin- 
djee,  1998,  Huber  and  Tsakmakis,  2000).  The  nonlinear  viscoelasticity  is  commonly  described 
through  a  rheological  spring-dashpot  model  that  features  a  rate-independent  equilibrium  part 
and  a  rate-dependent  viscous  part.  In  particular,  Zener  models  are  widely  used  (e.g.  Roland, 
1989,  Johnson  et  al.,  1995,  Bergstrom  and  Boyce,  1998,  Huber  and  Tsakmakis,  2000,  Quin- 
tavalla  and  Johnson,  2004,  Bergstrom  and  Hilbert,  2005,  Qi  and  Boyce,  2005,  Tomita  et  al., 
2008,  Areias  and  Matous,  2008).  Zener  models  are  composed  of  a  spring  (rate-independent 
part)  in  parallel  with  a  Maxwell  element  (rate-dependent  part).  The  constitutive  model  of  the 
time-dependent  part  includes  an  evolution  equation  for  the  viscous  deformation.  The  simplest 
evolution  equation  is  a  linear  viscous  flow  rule  (e.g.  Amin  et  al.,  2002,  Johlitz  et  al.,  2007). 
However,  experiments  have  shown  that  the  viscosity  at  finite  strains  is  a  function  of  the  driving 
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stress  (e.g.  Amin  et  al.,  2006),  the  strain  (e.g.  Amin  et  ah,  2006,  Hoo  Fatt  and  Ouyang,  2008) 
and/or  the  strain-rate  (e.g.  Khan  et  al.,  2006).  Thus,  nonlinear  viscous  flows  rules  have  been 
developed  using  a  micro-mechanism  inspired  approach  (e.g.  Bergstrom  and  Boyce,  1998,  Palm 
et  al.,  2006)  or  a  purely  phenomenological  approach  (e.g.  Khan  and  Zhang,  2001,  Colak,  2005, 
Amin  et  al.,  2006,  Khan  et  al.,  2006,  Hoo  Fatt  and  Ouyang,  2008). 

In  order  to  represent  the  equilibrium  contribution  as  the  nonlinear  elastic  spring,  rubber 
elasticity  models  are  employed  using  the  conventional  material  types  such  as  Neo-Hookean  (e.g. 
Johlitz  et  al.,  2007),  Mooney- Rivlin,  (e.g.  Huber  and  Tsakmakis,  2000),  or  Arruda  and  Boyce’s 
(1993)  eight-chain  model  (e.g.  Bergstrom  and  Boyce,  1998,  Tomitaet  al.,  2008).  Other  research 
groups  (e.g.  Amin  et  al.,  2002,  2006,  Hoo  Fatt  and  Ouyang,  2008)  used  modified  conventional 
material  models  to  describe  the  high  initial  stiffness  and  the  subsequent  strain  hardening. 

Although  many  researcher  have  recently  reported  on  experimental  behavior  of  various  poly¬ 
meric  materials  for  a  wide  range  of  strain  rates  (e.g.  Khan  and  Farrokh,  2006,  Mulliken  and 
Boyce,  2006,  Sarva  et  al.,  2007),  little  research  has  been  reported  on  the  modeling  of  elastomeric 
materials  for  loading  and  unloading  for  a  wide  range  of  strain  rates.  As  for  the  hereditary  inte¬ 
gral  approach  based  models,  most  multiplicative  decomposition  based  models  have  been  exper¬ 
imentally  validated  for  a  narrow  range  of  strain  rates.  Quintavalla  and  Johnson  (2004)  adopted 
the  Bergstrom-Boyce  model  to  describe  the  dynamic  behavior  of  cis-(l,4)  polybutadiene  at  high 
strain-rate  of  3000/s  to  5000/s.  Recently,  Hoo  Fatt  and  Ouyang  (2008)  proposed  a  thermo¬ 
dynamically  consistent  constitutive  model  with  modified  neo-Hookean  rubber  elastic  springs 
to  describe  the  steep  initial  stiffness  of  virgin  butadiene  rubber  under  tensile  and  compressive 
loading  at  intermediate  strain  rates  (76/s  to  450/s). 

The  present  work  focuses  on  the  modeling  of  the  loading  and  unloading  response  of  polyurea 
over  a  wide  range  of  strain  rates  (from  10~3/s  to  101/s)  and  at  large  compressive  strains  (up 
to  1.0).  Using  the  framework  of  the  multiplicative  decomposition  of  the  deformation  gradient, 
the  rheological  concept  of  two  parallel  Maxwell  elements  is  employed  to  develop  a  nonlinear 
viscoelastic  finite  strain  constitutive  model.  Guided  by  similar  developments  in  area  of  rate- 
dependent  plasticity  models  for  glassy  polymers  (e.g.  Qi  and  Boyce,  2005,  Anand  and  Ames, 
2006,  Anand  et  al.,  2009,  Ames  et  al.,  2009,  Ayoub  et  al.,  2009),  the  high  stiffness  of  elastomers 
at  small  deformations  is  modeled  using  Hencky’s  strain  energy  function.  In  view  of  the  simula- 
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tion  of  blast  and  impact  events,  we  limit  our  attention  to  the  mechanical  behavior  of  polyurea 
in  its  virgin  state.  After  presenting  all  experimental  results  in  Section  3.2,  the  constitutive 
model  is  detailed  in  Sections  3.3  and  3.4.  A  discussion  of  the  comparison  of  simulations  and 
experimental  results  for  polyurea  is  given  in  Section  3.5. 


3.2  Experimental  Investigation 

3.2.1  Material 

The  polyurea  specimens  used  in  the  study  are  extracted  from  a  12.7 mm  thick  polyurea  DragonShield- 
HT  Explosive  Resistant  Coating  (ERC)  layer  on  a  steel  armor  plate.  After  separation  from  the 
plates,  the  polyurea  specimens  of  diameter  Dq  =  10 mm  and  length  L$  =  10mm  are  machined 
using  conventional  milling  procedures.  The  shelf  life  of  the  polyurea  prior  to  testing  was  about 
4  years.  All  experiments  are  performed  on  polyurea  in  its  virgin  state  (condition  after  curing, 
no  prior  mechanical  loading  history). 

3.2.2  Relaxation  Experiments 

The  nonlinearity  of  the  viscoelastic  behavior  of  polyurea  at  low  and  intermediate  strain  rates 
is  identified  from  relaxation  tests.  Using  a  hydraulic  universal  testing  machine  (Model  8800, 
Instron),  the  material  is  rapidly  loaded  to  a  prescribed  strain  level  at  the  beginning  of  the 
experiment.  The  strain  is  subsequently  kept  constant  while  stress  relaxation  takes  place.  Nine 
relaxation  tests  are  performed  for  true  compressive  strains  ranging  from  So  —  0.048  to  1.0.  The 
engineering  strain  rate  during  the  rapid  loading  phase  was  about  0.48/s.  In  order  to  capture 
the  short-term  viscoelastic  behavior,  the  sampling  frequency  of  the  data  was  5 kHz;  the  stress 
is  recorded  for  about  10  minutes  to  identify  the  long-term  behavior.  Figure  3-la  shows  the 
recorded  stress  histories  for  various  constant  strain  loadings.  While  the  test  duration  was  not 
long  enough  to  observe  the  asymptotic  long-term  behavior  (i.e.  equilibrium  modulus),  the 
results  indicate  that  the  sampling  frequency  of  5kHz  is  high  enough  to  capture  the  asymptotic 
short-term  viscoelastic  behavior  (i.e.  initial  relaxation  modulus). 
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Figure  3-1:  Results  of  relaxation  tests  (left  column)  and  corresponding  simulations  (right  col¬ 
umn).  (a)-(b)  True  stress  histories,  (c)-(d)  Relaxation  moduli  histories,  (e)-(f)  Isochronous 
stress-strain  curves  for  different  instants  after  the  rapid  strain  loading. 
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The  relaxation  moduli  are  computed  using  the  definition 


=  —  (3.1) 

where  a  ( t )  is  the  time  history  of  the  measured  true  stress  for  a  constant  strain  of  So  (Figure 
3-la).  The  corresponding  relaxation  modulus  histories  (Figure  3-lc)  clearly  demonstrate  the 
nonlinear  viscoelastic  nature  of  the  material  response:  the  relaxation  modulus  decreases  from 
120 MPa  (at  the  strain  of  —0.048)  to  20 MPa  (at  the  strain  of  —1.0)  at  the  beginning  of 
relaxation;  after  more  than  10  minutes,  the  corresponding  relaxation  moduli  are  60 MPa  and 
10 MPa,  respectively.  The  strong  nonlinear  nature  of  the  viscoelastic  behavior  becomes  also 
apparent  when  plotting  the  isochronous  stress-strain  curves  (Figure  3-le). 

3.2.3  Continuous  Compression  Experiments 

In  addition  to  relaxation  experiments,  constant  strain  rate  loading  and  unloading  tests  are 
performed  at  absolute  true  strain  rates  between  1CT3  and  101/s  (Figure  3-2a).  The  measured 
true  stress  versus  logarithmic  strain  curves  for  five  different  strain  rates  are  summarized  in 
Figure  3-2b.  At  least  two  tests  are  performed  at  each  strain  rate  to  confirm  the  repeatability 
of  the  experimental  results.  During  the  loading  phase,  the  observed  stress  level  exhibits  strong 
strain  rate  sensitivity.  For  example,  the  observed  stress  level  for  10~3/s  is  about  twice  as  high 
as  for  101/s.  However,  during  unloading,  all  the  stress-strain  curves  seem  to  converge  to  the 
same  curve  regardless  of  the  unloading  strain  rate. 

3.2.4  Stair  Compression  Experiments 

In  order  to  investigate  the  stress  relaxation  behavior  during  loading  and  unloading,  experiments 
are  performed  with  stair-type  strain  history  profiles.  Figure  3-3a  shows  the  applied  strain  history 
as  a  function  of  the  normalized  time  t  |e/£o|-  The  same  maximum  strain  of  £o  =  —1-0  is  chosen 
for  all  experiments.  Experiments  are  performed  for  |eo|  =  1CT3/s,  10~2/s,  10-1/s,  10°/s  and 
101/s.  All  measured  stress-strain  curves  are  depicted  in  Figure  3-3b.  The  results  from  the  stair 
compression  tests  confirm  the  previous  observation  of  high  strain  rate  sensitivity  during  loading 
phase  along  with  nearly  rate  independent  behavior  during  unloading.  Also  the  magnitude  of 
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Figure  3-2:  Results  of  continuous  loading/unloading  experiments  and  corresponding  simula¬ 
tions.  (a)  History  of  applied  strain.  Note  that  time  in  the  x-axis  is  normalized  by  the  total 
testing  duration  |eo/£o|  where  £o  =  10-3  ,  10— 2 ,  10— 1 ,  10°,  and  10x/s  .  (b)  Stress-strain  curves 
obtained  from  experiments  and 
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stress  relaxation  during  the  loading  phase  is  found  to  be  larger  than  that  during  unlading.  The 
superposition  of  the  results  from  the  monotonic  and  stair  compression  experiments  (Figure  3- 
3d)  shows  that  the  specimens  recover  the  virgin  state  resistance  (described  by  the  curves  from 
the  monotonic  tests)  after  each  relaxation  period  in  the  stair  compression  experiment.  The 
resistance  recovery  appears  to  be  slower  for  large  strains  and  high  strain  rates. 

3.3  Constitutive  Model 

3.3.1  Motivation 

In  the  framework  of  multiplicative  decomposition  of  the  deformation  gradient,  models  of  the 
Zener  type  (Maxwell  element  in  parallel  with  a  spring)  are  widely  used  since  these  can  effec¬ 
tively  represent  the  nonlinear  viscoelastic  behaviors  of  many  elastomers.  For  Zener  models, 
it  is  common  to  identify  the  equilibrium  path  and  the  over-stress  contribution  from  the  ex¬ 
perimental  stress-strain  curves.  The  equilibrium  path  can  be  determined  in  an  approximate 
manner  from  two  alternative  approaches:  (1)  taking  the  mid-path  of  stress-strain  curves  from 
the  monotonic  loading  and  unloading  test  (e.g.  Bergstrom  and  Boyce,  1996),  (2)  connecting 
the  stress  relaxation  mid-points  between  loading  and  unloading  phase  from  the  stair-type  com¬ 
pression  tests  (e.g.  Qi  and  Boyce,  2005).  Figure  3-4  shows  the  equilibrium  paths  that  have 
been  identified  from  the  continuous  compression  tests  at  different  strain  rates.  The  summary 
plot  in  Figure  3-66a  clearly  shows  that  the  identified  equilibrium  paths  are  not  identical,  i.e. 
the  identified  equilibrium  path  depends  on  the  strain  rate.  For  instance,  at  a  strain  of  —1.0, 
an  equilibrium  stress  of  —10 MPa  is  found  from  the  experiment  at  10"3/s  while  an  equilibrium 
stress  of  —17 MPa  is  obtained  from  the  experiment  at  101/s.  This  important  observation  is 
confirmed  by  the  stair-type  compression  tests.  Recall  from  Figure  3-33b  that  the  amount  of 
stress  relaxation  is  different  for  loading  and  unloading.  Thus,  it  is  concluded  that  the  concept 
of  equilibrium  path  breaks  down  in  the  case  of  polyurea.  Consequently,  rheological  models  of 
the  Zener  type  will  not  be  able  to  describe  the  nonlinear  viscoelastic  behavior  of  polyurea  for 
both  loading  and  unloading. 

In  models  of  the  Zener  type,  the  equilibrium  part  is  represented  by  a  spring  element.  Qi  and 
Boyce  (2005)  explained  in  their  paper  on  the  micro-mechanism  inspired  constitutive  modeling 
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Figure  3-3:  Results  for  stair  compression  loading/unloading,  (a)  History  of  applied  strain. 
Note  that  time  in  the  x-axis  is  normalized  by  the  total  testing  duration  \£q/£q\  where  £q  — 
1  and  £o  =  10-3  ,  10 — 2 ,  10-1,  10°,  and  lO1/^  .  (b)  Stress-strain  curves  from  experiments 
and  (c)  Corresponding  simulations,  (d)  Comparison  of  results  from  continuous  and  stair-type 
loading/unloading  experiments  and  (e)  Corresponding  simulations. 
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Figure  3-4:  Estimation  of  the  behavior  of  Network  A  using  the  equilibrium  path  concept  for 
path  concept  for  different  strain  rates:  (a)  10_3/s,  (b)  10_2/s,  (c)  10_1/s,  (d)  10°/s  and  (e) 
lOVs 
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Figure  3-5:  Proposed  rheological  model  for  polyurea  composed  of  two  Maxwell  elements  in 
parallel. 

of  thermoplastic  polyurethanes  that  the  spring  element  associated  with  the  equilibrium  path 
represents  the  entropic  resistance  of  the  soft  domains  in  polyurethanes,  while  the  rate  dependent 
behavior  of  the  hard  part  is  represented  by  the  Maxwell  element  of  the  Zener  model.  Due  to 
the  apparent  rate  dependency  of  the  calculated  equilibrium  paths,  we  substitute  the  spring  of 
the  Zener  model  by  another  Maxwell  element.  The  resulting  rheological  model,  two  Maxwell 
springs  in  parallel  (Figure  3-5),  has  already  been  considered  in  the  past  by  Boyce  et  al.  (2000) 
and  Dupaix  and  Boyce  (2007)  to  describe  the  finite  strain  behavior  of  amorphous  polymers.  In 
the  context  of  amorphous  polymers  near  glass  transition  temperature,  two  Maxwell  elements 
have  been  used  to  represent  the  rate-dependent  inter-  and  intra-molecular  network  resistances 
under  monotonic  loading  conditions. 

In  the  following,  we  outline  the  constitutive  equations  for  each  of  the  four  basic  elements 
of  our  rheological  model.  Throughout  our  presentation,  we  refer  to  the  Maxwell  element  that 
is  predominantly  associated  with  the  deformation  resistance  of  the  soft  domain  of  polyurea 
as  “ Network  A”,  while  the  Maxwell  element  associated  with  the  effect  of  the  hard  domain  is 
referred  to  as  “ Network  B” . 

3.3.2  Homogenization 

The  macroscopic  deformation  gradient  Ftot  is  decomposed  into  a  rotation  free  volumetric  part 

1 

Froi  =  J1/31  with  J  =  det  Fj0t  (3.2) 
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and  an  isochoric  part  F 

F  =  J~l/3Ftot  (3.3) 

The  constitutive  model  for  the  stresses  induced  by  isochoric  deformation  is  developed  with 
the  Taylor  assumption  in  mind,  i.e.  the  isochoric  deformation  gradients  within  both  networks, 
Fyi  and  F b,  equal  the  macroscopic  isochoric  deformation  gradient 

F  =  FA=  F  b  (3.4) 


while  the  deviatoric  part  of  the  macroscopic  Cauchy  stress  T  corresponds  to  the  spatial  average 
of  the  local  stress  fields 


devT 


_1 

n 


/  crdflA  +  /  crdQ.s 

J  A  '  Jb 


(3.5) 


where  a  describes  the  local  deviatoric  stress  fields,  Q  is  the  total  volume  of  the  representa¬ 
tive  volume  element,  and  Q.A  and  Ob  are  the  corresponding  volumes  of  Networks  A  and  B, 
respectively.  In  the  context  of  micromechanism-inspired  polymer  models,  the  role  of  the  mi¬ 
crostructure  is  typically  neglected.  Instead,  the  constitutive  equations  for  the  individual  phases 
predict  a  weighted  macroscopic  Cauchy  stress  (in  the  present  case  the  weighted  deviator  stresses 
Ta  and  Tb)  such  that  the  macroscopic  stress  deviator  can  be  written  as 


devT  =  Ta  +  T  b 


(3.6) 


3.3.3  Constitutive  Equations  for  Volumetric  Deformation 

The  volume  change  of  polyurea  under  mechanical  loading  is  very  small  as  compared  to  the 
magnitude  of  isochoric  deformation.  For  simplicity,  we  assume  a  linear  elastic  relationship 
between  the  logarithmic  volumetric  strain  and  the  hydrostatic  part  of  the  macroscopic  Cauchy 
stress  tensor 

trT  In  J  , 

—  =  Krr  (3'7) 

where  k  denotes  the  bulk  modulus. 
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3.3.4  Maxwell  Model  for  Isochoric  Deformation 


The  responses  of  Network  A  and  Network  B  are  described  through  Maxwell  models.  In  this 
subsection,  we  provide  the  general  framework  for  the  formulation  of  the  constitutive  equations 
for  Maxwell  models.  In  the  Sections  3.3.5  and  3.3.6,  we  will  then  specialize  our  equations  for 
the  individual  networks. 

Kinematics 

Consider  a  Maxwell  element  K  subject  to  a  isochoric  deformation  gradient  F  k  with  det  F#  =  1. 
In  the  context  of  finite  strain,  we  assume  that  the  total  isochoric  deformation  gradient  can  be 
multiplicatively  decomposed  into  an  elastic  part  FeK  and  a  viscous  part  F^, 

F K  =  FrFk  (3-8) 

Furthermore,  we  introduce  the  total,  elastic  and  viscous  rate  of  isochoric  deformation  tensors 
Fk  =  LkFk  ,  FeK  =  LeKFeK  and  Ffc  =  LfcFfc  (3.9) 

with  the  relationship 

L*  =  14  +  F^LpV  (3.10) 

It  is  assumed  that  both  the  elastic  and  viscous  spins  are  zero  for  the  deformation  of  the  Maxwell 
elements.  Formally,  we  write 

Df,  :=  14  =  Lf  and  JDVK  :=  Lfc  =  L  f  (3.11) 

Due  to  the  spin- free  elastic  deformation,  the  rate  of  deformation  tensor  T)eK  may  also  be  ex¬ 
pressed  as  a  function  of  the  time  derivative  of  the  right  Cauchy-Green  tensor 

D  VK  =  \Y(KTCcKFeKl  (3.12) 

In  addition  to  the  tensor  description  of  the  kinematics,  we  make  use  of  the  scalar  deformation 
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measure 


(3.13) 


in  our  model  formulation.  The  strain-like  variable  £  is  zero  in  the  undeformed  configuration 
and  always  positive  for  deformed  configurations. 

Thermodynamics  and  Hyperelasticity 

For  isothermal  conditions,  the  second  law  of  thermodynamics  for  each  Maxwell  element  reads 

JTK  :  LK  —  ,ipK  >  0  (3.14) 

with  the  deviatoric  Cauchy  stress  T k  and  the  free  energy  ipK  (per  unit  initial  volume).  We 
impose  the  assumption  of  Green  elasticity  through  the  hyperelastic  relationship  (e.g.  Ogden, 
1984) 

T*  =  MF*^F'/}  (3‘15) 

where  CeK  =  F^ denotes  the  right  Cauchy-Green  tensor  of  isochoric  elastic  deformation. 
Recall  that  we  assume  isochoric  elastic  deformation  which  implies  detF^  =  1  and  trD^  =  0. 
Using  Eq.  (3.10)  and  Eq.  (3.15)  in  Eq.  (3.14)  yields  the  thermodynamic  constraint 

JTk  ■  (F^L^F^1)  =  (jFf  T*F^t)  :  Dfc  >  0  (3.16) 

In  the  following,  the  Mandel  stress 

M K  :=  Jdev  {f£tkF f}  (3.17) 

is  used  as  the  driving  stress  of  viscous  deformation  while  the  viscoplastic  constitutive  equations 
are  written  such  that 

Mk  :Dvk>0  (3.18) 
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Flow  Rule  for  Viscous  Flow 


The  direction  of  viscous  flow  is  supposed  to  be  aligned  with  the  direction  of  the  driving  Mandel 
stress.  Formally,  we  write  the  flow  rule 


Ma'  =  VK&K 


(3.19) 


with  the  deformation  and  rate  dependent  viscosity  rjK  >  0.  The  rate  dependent  viscosity  is 
defined  through  the  nonlinear  relationship  rhk  =  rnjt  (<4)  between  the  equivalent  Mandel  stress 


mk 


-Mfc  ; 


(3.20) 


and  the  equivalent  viscous  strain  rate 

dk  =  :  D - 


(3.21) 


Using  Eq.  (3.20)  and  Eq.  (3.21)  in  Eq.  (3.19),  we  obtain 


Vk 


2  mk 

3  dk 


(3.22) 


3.3.5  Specialization  of  the  Maxwell  Model  for  Network  A 

Network  A  represents  the  rubbery  response  of  the  soft  domain  which  is  modeled  through  a 
Maxwell  element  composed  of  a  nonlinear  elastic  Gent  spring  and  a  nonlinear  viscous  damper 
that  accounts  for  the  apparent  rate  dependency  of  the  theoretical  equilibrium  path. 


Gent  Elasticity  of  Network  A 

For  isochoric  deformation,  Gent’s  (1996)  free  energy  function  may  be  written  as 

ip  A  =  -\vaJa  In  ^1  ~  tTCjA  3)  (3-23) 
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with  the  material  parameters  ^A>  0  (initial  modulus)  and  J/ 1  >  0  (locking  stretch).  With  this 
form  of  free  energy,  the  Cauchy  stress  TA  in  Network  A  reads: 

trC}A~3)  devB eA  (3.24) 

where  BA  —  denotes  the  right  Cauchy-Green  tensor  of  isochoric  elastic  deformation. 

Nonlinear  Viscous  Response  of  Network  A 

Figure  3-6a  suggests  that  the  viscoelastic  contribution  to  the  stress  of  Network  A  increases 
monotonically  as  the  deformation  increases.  Furthermore,  the  experimental  results  suggest  a 
power-law  relation  between  the  viscosity  rjA  and  the  equivalent  rate  of  viscous  deformation. 
Thus,  the  nonlinear  viscous  evolution  law  for  Network  A  is  written  as 

mA  =  ^r]AdA  =  PA{exp(-l)  (3'25) 

with  the  reference  rate  of  deformation  dref  =  1/s,  the  viscosity  constant  PA  >  0,  and  the 
exponent  nA>  0  controlling  the  rate-sensitivity. 

3.3.6  Specialization  of  the  Constitutive  Equations  for  Network  B 

Network  B  represents  the  response  of  the  hard  domain  which  is  modeled  through  a  Maxwell 
element  composed  of  a  linear  elastic  spring  and  a  nonlinear  viscous  damper  that  accounts  for 
the  high  stiffness  at  small  strain  and  the  rate  sensitivity  of  hysteresis. 

Hencky  Elasticity  of  Network  B 

The  stiffness  of  Network  B  governs  the  stiff  macroscopic  response  of  polyurea  at  small  strains. 
Here,  we  make  use  of  the  isochoric  part  of  Hencky’s  strain  energy  function 

V'b  =  ^/Ab  (InU^j)  :  (lnUfj)  (3.26) 

where  the  elastic  stretch  tensor  \JeB  is  found  from  the  polar  decomposition  of  the  elastic  defor¬ 
mation  gradient,  FeB  =  R^U#  with  R^R#  =  1.  Upon  evaluation,  we  obtain  the  deviatoric 
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part  of  the  Cauchy  stress  acting  on  Network  B 


Tb 


J 


(lnU|)Rf 


(3.27) 


Nonlinear  Viscous  Response  of  Network  B 


Figure  3-7a  suggests  that  the  overall  viscoelastic  response  of  Network  B  is  similar  to  the  linear 
viscoelastic  response  of  a  Maxwell  model  composed  of  a  linear  spring  (of  Young’s  modulus  E) 
and  a  linear  viscous  damper  (of  viscosity  rj).  The  stress-strain  (cr  —  s)  response  of  a  Maxwell 
model  with  the  linear  components  at  a  constant  strain  rate  reads  (e.g.  Brinson  and  Brinson, 


2008) 


a  —  Eteq 


1  —  exp 


(3.28) 


where  r  =  rj/E  is  the  relaxation  time,  and  Eq  is  the  applied  constant  strain  rate.  Note  that 
even  though  the  Maxwell  model  is  composed  of  only  linear  elements,  the  stress-strain  response 
is  nonlinear  at  constant  strain  rates.  However,  our  experiments  suggest  that  the  stress-strain 
response  of  Network  B  is  linear  for  small  strains  regardless  of  strain  rate. 

Inspired  by  the  results  for  the  one-dimensional  Maxwell  model  with  linear  components 
(Eq.  (3.28)),  we  propose  a  constitutive  equation  that  guarantees  a  linear  elastic  response  of 
Network  B  at  small  strains 


mB 


VShbQb 


d-B 

dref 


tib 


1  —  exp 


Qb 


(3.29) 


where  the  dimensionless  material  constants  Qb  and  ns  characterize  the  viscous  properties  of 
Network  B. 


3.4  Identification  of  the  Model  Parameters 

The  proposed  constitutive  model  comprises  eight  material  parameters:  four  parameters  (fxA,  Ja, 
Pa,  ua)  for  Network  A,  three  parameters  (p, B,  Qb,  ns)  for  Network  B ,  and  one  parameter  (k) 
describing  the  elastic  volumetric  response.  The  constitutive  equations  are  detailed  for  uniaxial 
stress  loading  in  order  to  identify  all  material  parameters  from  the  monotonic  loading  and 
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unloading  compression  experiments. 

3.4.1  Constitutive  Equations  for  Uniaxial  Loading 

Under  uniaxial  stress  loading,  the  macroscopic  stress  tensor  takes  the  form 

T  =  ae i  ®  ei  (3.30) 


where  a  is  the  true  stress  (macroscopic  Cauchy  stress)  while  the  macroscopic  deformation 
gradient  reads 

Ff at  =  Aiei  ®  ex  +  A2  (e2  ®  e2  +  e3  <8>  e3)  (3.31) 

J  =  detFtot  =  \\  (A2)2  (3.32) 


At  the  network  level,  the  isochoric  part  of  the  deformation  gradient  may  be  expressed  as  a 
function  of  the  stretch  A, 


F  =  Aei  ®  ei  +  —j=  (e2  ®  e2  +  e3  <g>  e3) 
v  A 


with  A  = 


(3.33) 


The  total  stretch  of  a  Maxwell  element  K  can  be  decomposed  further  into  its  elastic  part  XeK 
and  the  viscous  part  XVK 

A  =  XeKXvK  (3.34) 


Using  this  notation,  the  kinematic  variables  used  in  our  constitutive  equations  may  be 
expressed  as  a  function  of  the  elastic  and  viscous  stretches  and  their  respective  time  derivatives 


trC 


e  _ 

K  ~ 


(3.35) 


devB^  -  l 


(Afc)2  -  Te~ 


K  J 


(2ei  ®  ei  -  e2  <8»  e2  -  e3  ®  e3) 


InU eK  =  -  (In  XeK)  (2ei  ®  ex  -  e2  ®  e2  -  e3  ®  e3) 


(3.36) 

(3.37) 


C  = 


(3.38) 
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(3.39) 


The  deviatoric  stress  tensors  T k  acting  on  the  Maxwell  elements  are  written  as 


T k  =  -!r  (2ei  ®  ei  -  e2  ®  e2  -  e3  (g)  e3) 
o 


(3.40) 


C  —  O'  A  +  &  B 


(3.41) 


a  a  and  oB  are  true  stresses  acting  on  Networks  A  and  B,  respectively. 

Using  expressions  Eq.  (3.30)  to  Eq.  (3,41),  we  obtain  the  simplified  constitutive  equations 
for  uniaxial  loading: 

•  Volumetric  deformation  induced  linear  elasticity  with  the  model  parameter  k 


a  —  3k- 


(3.42) 


Gent  elasticity  of  Network  A  with  the  model  parameters  Ja  and 


AM  f/xe  x 2  1  1  1  f/xe  x2  ,  2 


aA=j  TsFM-1 


(3.43) 


Nonlinear  viscous  response  of  Network  A  with  the  model  parameters  Pa  and  ua 


M  =  ^(expC-l)  (AAj 


(3.44) 


Hencky  elasticity  of  Network  B  with  the  model  parameter  /iB 


Nonlinear  viscous  response  of  Network  B  with  the  model  parameters  and  , 


(3.45) 
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Q  B  \  dref 


(3.46) 
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Table  3.1:  Summary  of  material  parameters  indentified  from  monotonic  loading/unloading  tests 
under  five  different  strain  rates. _ 


Rubber 

Pa 

7.00  MPa 

Isochoric  Part 

Spring 

Ja 

10.7 

of  Network  A 

Viscous 

Pa 

4A2MPa 

Damper 

Ua 

0.0646 

Linear 

Spring 

Jb 

82.3  MPa 

Isochoric  Part 
of  Network  B 

Viscous 

Qa 

0.0447 

Damper 

nA 

0.0755 

Volumetric  Part 

K 

829 MPa 

3.4.2  Model  Calibration 

The  1-D  calibration  process  is  composed  of  two  steps:  the  seven  parameters  (i.e.  /.iA,  Ja, 
P/ 1,  nA ,  [iB,  Qb  and  nB)  relating  to  isochoric  deformation  are  identified  assuming  material 
incompressibility  (J  =  1),  before  the  bulk  modulus  k  is  estimated  based  on  the  Lame  moduli 
and  the  Poison’s  ratio.  A  summary  of  all  identified  material  parameters  is  given  in  Table  1. 

Material  Parameters  of  Network  A  (nA,  Ja ,  Pai  n A ) 

The  model  is  calibrated  such  that  the  response  of  Network  A  describes  the  rate  dependency 
of  the  hypothetical  equilibrium  path  in  the  continuous  compression  experiments.  As  discussed 
in  Section  3.3.1  (see  also  Figure  3-6a),  the  hypothetical  equilibrium  path  a  a  —  or  a  (A,  Ja)  for 
a  given  strain  rate  amplitude  dA  is  approximated  by  the  average  of  the  measured  loading  and 
unloading  stress-strain  curves  at  a  total  strain  rate  e 

a  a  (A,  dA)  ~  i  [a  A  (A,  dA )  +  crA  (A,  - dA )]  with  dA  ^  |e|  (3.47) 

The  material  parameters  for  Network  A  are  calibrated  in  an  iterative  manner: 
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1.  Assuming  A  A(t)  cs  A  (t),  we  estimate  (/j,a,JA)  from  the  best  fit  of  Eq.  (3.4344)  to  the 
experimentally  obtained  stress-strain  curve  for  the  highest  strain  rate  case  (here,  |e|  ~ 

lOVs). 

2.  For  the  remaining  strain  rates,  we  calculate  \eA  ( t )  based  on  Eq.  (3.43)  using  ()iA,  Ja)- 

3.  Using  the  relationship  A^  it)  =  A  ( t )  /XA  ( t ),  we  calculate  Icr^l  as  a  function  of  the  strain¬ 
like  variable  (  for  each  experiment  (Figure  3-6b);  in  the  figure,  the  curve  labels  indicate 
the  average  effective  viscous  strain  rate  dA . 

4.  The  parameters  PA  and  ha  are  then  obtained  from  the  approximation  Eq.  (3.44)  of  the 
curves  shown  in  Figure  3-6c. 

5.  It  is  useful  to  evaluate  the  term  IU  :=  \<ja\/  (expC  —  1)  from  each  experimental  curve. 
11,4  is  a  function  of  the  effective  viscous  strain  rate  dA  only  and  independent  of  the  total 
strain.  In  other  words,  each  curve  shown  in  Figure  3-6b  reduces  to  a  single  data  point 
in  Figure  3-6c.  The  subsequent  curve  fit  11^  =  (Pa/ J)  ( dA/dref)nA  of  the  data  points  in 
Figure  3-6c  yields  the  model  parameters  Pa  and  n^. 


After  completing  Steps  1  to  5,  we  obtain: 

HA  =  6.85 MPa  ,  JA  =  11.0 
PA  =  4.44MPa  ,  nA  =  0.0648 


(3.48) 


A  second  calibration  loop  is  performed,  starting  with  A^  (t)  =  A  (t) /A vA(t)  in  Step  1,  where 
A^  (t)  is  the  viscous  deformation  predicted  by  the  model  parameters  Eq.  (3.48).  Subsequent 
evaluation  yields: 


fiA  =  7.00 MPa  ,  Ja  =  10.7 
PA  =  4.42 MPa  ,  nA  =  0.0646 


(3.49) 


The  calibration  procedure  has  been  stopped  at  this  point  as  additional  iterations  did  not  improve 
the  accuracy  of  the  curve  fit. 
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Approximate  Network  A  Contribution 


Network  A:  Calibration  of  n. 
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Figure  3-6:  Material  model  parameter  identification  for  Network  A.  (a)  Estimated  experimen¬ 
tal  stress-strain  curve  for  Network  A\  (b)  Equivalent  Mandel  stress  as  a  function  of  the  scalar 
deformation  measure  (c)  II 4  as  function  of  the  equivalent  viscous  deformation  rate  dA  for 
the  identification  of  the  model  parameters  Pa  and  tia  through  power-law  fit. 
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Material  Parameters  of  Network  B  (frB,  QB,  ng) 

The  stress  contribution  of  Network  B  is  obtained  by  subtracting  the  identified  stress  con¬ 
tribution  of  Network  A  (using  the  parameters  Eq.  (3.49)  in  Eq.  (3.44))  from  the  original 
stress-strain  curve  for  each  strain  rate 

crB  (A,e)  =  a(A,e)  -  aA  (X,3a)  (3.50) 

where  3a  corresponds  to  the  average  viscous  strain  rates  that  have  been  determined  throughout 
the  calibration  of  Network  A.  Figure  3-7  shows  the  corresponding  results  for  different  strain 
rates.  According  the  proposed  rheological  model,  viscoelastic  behavior  of  Network  B  should 
be  symmetric  with  respect  to  the  stress,  i.e.  (A,  — e)  =  —  crB(X,  e).  However,  the  curves 

in  Figure  3-7  illustrate  that  this  symmetry  assumption  holds  only  approximately  true  when 
using  the  calibrated  analytical  expression  Eq.  (3.44)  for  Network  A.  In  close  analogy  with 
the  calibration  procedure  for  Network  A,  we  identify  the  model  parameters  associated  with 
Network  B: 

1.  We  assume  \eB  (t)  ~  A  (t)  for  the  highest  strain  rate  case  (i.e.  strain  rate  of  101/s). 

2.  Determine  jiB  from  the  best  curve  fit  for  the  experimental  stress-strain  curve  for  the 
loading  phase  using  Eq.  (3.45). 

3.  Using  fiB>  we  can  then  calculate  X%  ( t )  and  XVB  ( t )  =  A  (t)  /XeB  ( t )  for  the  other  experiments. 
The  corresponding  \aB\  versus  C  curves  are  shown  in  Figure  3- 7b.  According  to  Eq.  (3.46), 
we  have 

M  =  %B  1~exp(-^)  (3-51) 

where  (B  :=  QB  (dB/ dref)nB  is  a  function  of  the  strain  rate  only.  Thus,  we  determine 
(B  from  the  approximation  of  Eq.  (3.51)  to  the  |<7jg|  vs.  (  curves. 

4.  Subsequently,  the  material  constants  QB  and  nB  are  found  from  plotting  (Basa  function 
of  the  effective  viscous  strain  rate  dB. 
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A  second  parameter  identification  loop  is  performed  using  A^  (t)  =  A  (t)  /XVB  (t)  leading  to: 


=  82.3MPa 

Qb  =  0.0447  ,  nA  =  0.0755 


(3.52) 


Bulk  Modulus  of  Material  (k) 

The  initial  shear  modulus  for  the  proposed  model  is  the  sum  of  the  shear  modulus  of  each 
network,  i.e.,  fi  =  fiA  +  \iB  =  89. 3M Pa.  Experiments  revealed  that  the  Poisson’s  ratio  of 
v  —  0.45  is  nearly  rate- independent  throughout  the  deformation  up  to  the  strain  of  1.0  (Shim 
and  Mohr,  2009a).  Since  the  free  energy  for  the  volumetric  deformation  is  based  on  the  linear 
elastic  model,  we  can  determine  the  bulk  modulus  using  the  relation  between  the  elastic  moduli 

K=W^  =  S2mPa  (3'53) 

3.5  Comparison  of  Simulation  and  Experiments 

The  constitutive  model  has  been  programmed  as  a  user  material  subroutine  for  the  finite  element 
software  ABAQUS/Explicit.  Using  the  identified  material  model  parameters,  we  perform  the 
simulations  of  all  experiments  described  in  Section  3.2. 

3.5.1  Continuous  Compression 

Figure  3-2c  summarizes  the  numerically  predicted  stress-strain  curves  for  different  strain  rates. 
Similarly  to  the  experiments  (Figure  3- 2b),  the  simulation  results  exhibit  a  high  strain  rate 
sensitivity  during  the  loading  phases.  In  addition,  the  characteristic  convergence  of  the  stress- 
stain  curves  for  different  loading  velocities  is  also  captured  by  the  simulations.  The  comparisons 
of  experiment  and  simulation  for  individual  strain  rates  are  shown  in  Figures  3-8a  to  3-8e.  The 
stress-strain  curves  show  a  better  agreement  during  loading  than  during  unloading.  The  loading 
portions  typically  feature  the  largest  error  for  small  strains  (less  than  0.2)  while  good  predictions 
are  achieved  for  large  strains.  The  predicted  response  during  unloading  is  too  stiff  for  all  strain 
rates,  i.e.  the  simulations  reach  a  stress  of  zero  at  larger  compressive  strains  than  measured  in 
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Approximate  Network  B  Contribution 


Network  B:  Calibration  of  nB  &  GB 
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Figure  3-7:  Material  model  parameter  identification  for  Network  B.  (a)  Estimated  experi¬ 
mental  stress-strain  curve  for  Network  B ;  (b)  Equivalent  Mandel  stress  as  a  function  of  the 
scalar  deformation  measure  (c)  £ B  as  function  of  the  equivalent  viscous  deformation  rate  rig; 
identification  of  the  model  parameters  Qb  and  rig- 
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the  experiments.  From  an  energetic  point  of  view,  it  is  worth  noting  that  the  differences  during 
the  loading  phase  underestimate  the  mechanical  work,  while  the  work  is  overestimated  during 
unloading.  However,  these  two  errors  compensate  each  other  when  calculating  the  viscous 
energy  dissipation  for  the  entire  loading/unloading  cycle.  In  other  words,  despite  small  errors 
in  the  stress  level  during  loading  and  unloading,  the  overall  energy  dissipation  is  well  predicted 
over  a  wide  range  of  strain  rates. 

The  graphs  in  Figure  3-8  also  depict  the  predicted  individual  contributions  of  Network  A 
(dotted  blue  line)  and  Network  B  (dashed  green  line).  During  the  model  calibration,  we 
assumed  that  the  effective  viscous  strain  rate  dj\  is  the  same  during  loading  and  unloading, 
which  implies  that  Network  A  exhibits  the  same  stress-strain  response  during  loading  and 
unloading.  However,  since  the  total  viscous  strain  keeps  on  decreasing  during  unloading  (the 
driving  stress  for  Network  A  is  compression  during  both  loading  and  unloading)  even  though 
A  increases,  the  decrease  of  the  stress  level  is  faster  during  unloading  than  the  corresponding 
increase  during  loading  (hysteresis  behavior  of  a  Maxwell  element).  The  differences  between 
the  loading  and  unloading  portion  for  Network  A  (blue  dashed  curves)  illustrate  this  behavior. 

The  model’s  symmetric  response  of  Network  B  (the  green  dashed  curves  show  only  the 
loading  portion)  is  characterized  by  a  plateau  over  a  wide  range  of  strains.  Since  the  elastic 
strains  in  Network  B  are  small,  the  driving  stress  of  viscous  deformation  in  Network  B  changes 
rapidly  from  compression  during  loading  to  tension  during  unloading.  As  a  result,  the  effective 
rate  of  viscous  deformation  cIb  is  the  same  during  loading  and  unloading  which  is  consistent  with 
the  assumptions  made  during  calibration.  Thus,  it  is  concluded  that  the  differences  between 
the  experiments  and  model  predictions  are  associated  with  the  response  of  Network  A. 

3.5.2  Stair  Compression 

Figure  3-3c  provides  an  overview  on  the  simulation  results  for  all  loading  cases  next  to  the 
experimentally-measured  stress-strain  curves  from  stair  compression  tests.  A  direct  comparison 
of  the  experimental  and  the  simulation  results  for  each  loading  velocity  is  shown  in  Figure  3- 
9a  to  3-9e.  Overall,  we  observe  the  same  model  accuracy  as  for  the  monotonic  loading  and 
unloading  experiments,  i.e.  the  model  slightly  underestimates  the  stress  level  for  small  strains 
and  predicts  a  faster  decrease  in  stress  level  after  reversal  of  the  loading  direction. 
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Figure  3-8:  Comparison  of  simulation  results  and  experiments  for  continuous  loading-unloading 
cycles,  (a)  10_3/s,  (b)  10_2/s,  (c)  10-1/s>  (d)  10°/s  and  (e)  lOVs 
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Figure  3-9:  Comparison  of  simulation  results  and  experiments  for  stain  loading/unloading,  (a) 
10_3/s,  (b)  10_2/s,  (c)  10-1/s>  (d)  10°/s  and  (e)  lO^s 
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The  relaxation  behavior  of  Network  B  is  symmetric  for  loading  and  unloading  (green  dashed 
lines  in  Figure  3-9).  Conversely,  the  stress-strain  curves  for  Network  A  (blue  dotted  lines  in 
Figure  3-9)  reveal  that  stress  relaxation  is  more  pronounced  during  the  loading  phase,  while 
unloading  is  mainly  elastic.  This  non-symmetric  feature  of  Network  A  is  also  visible  at  the 
macroscopic  level  where  the  model  predicts  the  experimental  observation  of  a  faster  relaxation 
during  loading  than  during  unloading.  The  model  provides  accurate  estimates  of  the  relaxation- 
induced  stress  decrease  (jump)  for  small  strains  during  the  loading  phase,  while  this  decrease 
is  underestimated  at  large  strains.  This  is  attributed  to  the  fact  that  the  stress  contribution 
of  Network  B  reaches  its  plateau  level  at  small  strains;  as  a  result,  the  amplitude  of  stress 
relaxation  within  Network  B  is  very  similar  for  small  and  large  strains. 

3.5.3  Relaxation 

Figure  3- lb  shows  the  predicted  stress  histories  for  all  relaxation  experiments.  In  addition, 
the  histories  of  the  relaxation  modulus  are  presented  in  Figure  3- Id.  The  relaxation  modulus 
histories  are  well  predicted  for  the  entire  range  of  strains.  However,  the  comparison  of  the 
true  stress  versus  time  curves  indicates  that  the  model  predicts  slower  stress  relaxation  than 
observed  in  the  experiments.  Similar  conclusions  may  be  drawn  from  the  comparison  of  the 
experimental  and  simulated  isochronous  stress- strain  curves.  Observe  that  the  instantaneous 
stress-strain  curves  agree  well  with  the  test  results,  while  the  long-term  behavior  (e.g.  blue 
curves  for  616s)  predicted  by  the  model  overestimates  the  stress  level  at  large  strains. 

3.5.4  Discussion 

The  model  is  able  to  predict  the  strain  rate  dependent  loading  and  relaxation  behavior  polyurea 
with  reasonable  accuracy.  Furthermore,  the  model  provides  good  estimates  of  the  rate  depen¬ 
dent  viscous  dissipation  throughout  the  loading/unloading  cycles.  In  the  present  rheological 
framework  of  two  parallel  Maxwell  elements,  the  identification  of  the  behavior  of  Network  A 
is  important.  Here,  the  response  of  Network  A  has  been  defined  through  the  “rate-dependent 
equilibrium  path”  concept;  the  mechanical  response  of  Network  A  is  identical  throughout  com¬ 
pressive  loading  and  compressive  unloading.  Furthermore,  in  the  particular  case  of  polyurea, 
this  imaginary  equilibrium  path  depends  on  the  absolute  value  of  the  total  true  strain  rate.  The 
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response  of  the  calibrated  constitutive  model  illustrates  that  the  assumption  of  a  Maxwell  model 
does  not  lead  to  an  exact  description  of  the  material  behavior  for  both  loading  and  unloading. 
The  hysteresis  behavior  for  loading  and  unloading  at  a  constant  total  strain  rate  is  an  inherent 
feature  of  Maxwell  models.  Consequently,  the  rate  dependency  of  the  conceptual  equilibrium 
path  (which  is  hysteresis  free)  cannot  be  represented  accurately  by  a  Maxwell  model.  An  at¬ 
tempt  has  been  made  to  calibrate  the  material  model  parameters  using  a  numerical  parameter 
identification  scheme  (minimization  of  the  error  between  the  experiments  and  simulations) .  In 
other  words,  the  equilibrium  path  concept  based  separation  of  the  contributions  of  Network  A 
and  Network  B  has  been  omitted  throughout  the  material  model  parameters  identification. 
However,  the  identified  model  response  was  not  more  satisfactory  as  improvements  on  the  un¬ 
loading  response  could  only  be  obtained  at  the  expense  of  the  quality  of  the  loading  curve 
prediction. 

In  order  to  come  up  with  an  improved  constitutive  model,  one  may  consider  a  different  rheo¬ 
logical  model  for  Network  A,  but  to  the  best  of  the  authors’  knowledge,  no  thermodynamically 
consistent  model  has  been  published  in  the  open  literature  that  can  model  the  apparent  elastic 
(same  loading  and  unloading  curve),  but  yet  strain  rate  dependent  behavior  of  Network  A. 
Thus,  alternative  modeling  approaches  need  to  be  explored  in  the  future.  One  thermody¬ 
namically  consistent  option  is  the  introduction  of  loading  and  unloading  conditions  for  viscous 
evolution.  However,  more  research  is  needed  at  the  microstructural  (i.e.  molecular)  level  to 
justify  such  loading  and  unloading  conditions  from  a  mechanism  point  of  view. 

3.6  Conclusions 

The  large  strain  compression  response  of  polyurea  is  investigated  for  strain  rates  ranging  10_3/s 
from  to  101  / s.  Continuous  and  stair-like  loading  and  unloading  experiments  are  performed  in 
addition  to  simple  relaxation  tests.  The  experimental  results  reveal  that  the  so-called  equilib¬ 
rium  path  concept  breaks  down  in  the  case  of  polyurea.  In  contradiction  with  the  definition  of 
the  equilibrium  path,  the  identified  average  stress-strain  curves  in  continuous  loading/unloading 
experiments  depend  on  the  rate  of  loading.  Thus,  as  an  alternative  to  the  equilibrium  path 
based  rheological  models  of  the  Zener  type,  a  new  constitutive  model  is  formulated  assuming 
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two  parallel  Maxwell  elements.  The  finite  strain  constitutive  equations  are  outlined  in  detail. 
Subsequently,  the  eight  material  model  constants  are  calibrated  to  describe  the  mechanical  be¬ 
havior  of  polyurea.  The  model  predictions  are  in  good  agreement  with  the  experimental  results. 
The  characteristic  convergence  of  the  unloading  paths  for  different  strain  rates  is  successfully 
captured  by  the  model.  Furthermore,  the  rate  dependent  viscous  dissipation  throughout  a  load¬ 
ing/unloading  cycle  is  described  with  reasonable  accuracy.  The  predicted  unloading  behavior 
is  too  stiff  which  is  associated  with  the  limitations  of  the  Maxwell  model  that  describes  the  soft 
domain  of  polyurea. 
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Chapter  4 


Validation  Application 

4.1  Introduction 

Polyurea  is  a  highly  viscoelastic  rubber  material  that  is  used  for  the  impact  protection  of  vehicle 
structures.  It  is  considered  for  the  armor  protection  and  retrofitting  of  military  vehicles  that  are 
exposed  to  the  blast  loading  of  improvised  explosive  devices.  The  anticipated  effect  of  polyurea 
coatings  on  the  blast  resistance  of  steel  plates  is  twofold.  Firstly,  the  polyurea  can  directly 
absorb  a  portion  of  the  blast  energy  as  it  undergoes  large  deformations.  Secondly,  the  onset  of 
ductile  fracture  of  a  steel  plate  may  be  retarded  through  the  use  of  a  polyurea  coating,  thereby 
increasing  the  energy  absorption  of  the  steel  structure. 

As  discussed  by  Xue  and  Hutchinson  (2008),  necking  occurs  under  uniaxial  tension  when 
the  average  true  stress  becomes  equal  to  the  overall  tangent  hardening  modulus  (Considere 
criterion).  In  the  case  of  a  coated  ductile  substrate,  a  high  strain  hardening  coating  material 
can  increase  the  effective  hardening  modulus  of  the  bilayer  material  such  that  necking  is  re¬ 
tarded  with  respect  to  the  Considere  strain  of  the  uncoated  material.  The  bifurcation  analysis 
of  Guduru  et  al.  (2006)  reveals  that  an  added  surface  layer  can  increase  the  resistance  of  a 
structural  element  to  fragmentation.  Moreover,  their  results  show  that  the  addition  of  a  soft 
coating  with  high  strain  hardening  can  improve  the  weight  specific  energy  absorption  of  the 
structural  element.  Xue  and  Hutchinson  (2008)  demonstrate  that  the  ratio  of  the  elastomer 
modulus  to  the  flow  strength  of  the  substrate  controls  the  effect  of  necking  retardation.  Mc- 
Shane  et  al.  (2008)  performed  tension  and  bulge  tests  on  copper/polyurethane  bilayers  under 
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static  and  dynamic  conditions.  Their  experimental  measurements  indicate  that  coatings  do  not 
provide  dynamic  performance  benefits  on  an  equal  mass  basis.  While  the  total  blast  resistance 
increases,  the  weight  specific  energy  absorption  of  the  structure  may  actually  decrease  through 
the  application  of  a  polymer  coating.  Dynamic  ring  expansion  experiments  have  been  performed 
by  Zhang  et  al.  (2009)  on  polyurea  coated  aluminum  6061  —  O  and  copper  101  at  very  high 
strain  rates  (4000  —  15000/s).  Their  experimental  results  show  that  there  is  no  significant  effect 
of  the  polyurea  coating  on  the  strain  at  the  onset  of  localization. 

It  appears  that  the  neck  retardation  effect  in  coated  ductile  substrates  is  difficult  to  achieve 
when  using  polyurea  in  combination  with  typical  engineering  materials.  However,  as  pointed 
out  by  McShane  et  al.  (2008),  polyurea  coatings  may  still  be  seen  as  a  practical  solution  for 
enhancing  the  blast  resistance  of  metallic  structures  because  of  the  ease  of  applying  polyurea 
on  existing  structures  (retrofitting).  Even  though  the  performance  of  the  steel  substrate  may 
remain  unaffected,  a  very  thick  poly  urea  layer  can  still  increase  the  energy  absorption  in  absolute 
terms.  The  impulsive  loading  experiments  of  Amini  et  al.  (2009a)  reveal  that  polyurea  coatings 
have  a  strong  effect  on  the  energy  transfer  to  the  steel  plate.  In  particular,  they  demonstrate 
that  the  positioning  of  polyurea  on  the  impact  side  promotes  failure  of  the  steel  plate  under 
shock  loading  while  a  polyurea  layer  on  the  back  of  the  plate  attenuates  the  shock.  In  the  present 
paper,  we  deal  with  the  prediction  of  the  large  deformation  behavior  of  polyurea  in  structural 
applications.  Xue  and  Hutchinson  (2008)  made  use  of  a  Mooney-Rivlin  model  for  the  polymer 
coating  in  their  numerical  analysis  of  the  polymer/metal  bilayers.  Zhang  et  al.  (2009)  modeled 
the  behavior  of  polyurea  using  a  nonlinear  hyperelastic  material  model.  However,  both  uniaxial 
compression  and  tension  tests  have  demonstrated  that  the  mechanical  response  of  polyurea  is 
highly  strain-rate  dependent  (e.g.  Amirkhizi  et  al.,  2006,  Roland  et  al.,  2007,  Sarva  et  al.,  2007, 
Shim  and  Mohr,  2009a).  Amini  et  al.  (2009b)  make  use  of  the  temperature-,  rate-  and  pressure- 
sensitive  constitutive  model  by  Amirkhizi  et  al.  (2006)  to  provide  supporting  simulation  results 
of  their  direct  pressure  pulse  experiments. 

Finite  viscoelasticity  models  of  elastomers  may  be  formulated  using  the  so-called  hereditary 
integral  approach  (Coleman  and  Noll,  1961,  Bernstein  et  al.,  1963,  Lianis,  1963,  McGuirt  and 
Lianis,  1970,  Leonov,  1976,  Johnson  et  al.,  1994,  Haupt  and  Lion,  2002,  Amirkhizi  et  al.,  2006) 
but  their  validity  is  often  limited  to  a  narrow  range  of  strain  rates  (Yang  et  al.,  2000,  Shim  et  al., 
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2004,  Hoo  Fatt  and  Ouyang,  2007).  As  an  alternative  to  the  hereditary  integral  approach,  the 
framework  of  multiplicative  decomposition  of  the  deformation  gradient  (Kroner,  1960  and  Lee, 
1969)  is  frequently  used  in  finite  viscoelasticity  (e.g.  Sidoroff,  1974,  Lubliner,  1985,  Le  Tallec 
et  ah,  1993,  Reese  and  Govindjee,  1998,  Huber  and  Tsakmakis,  2000).  In  that  framework, 
the  nonlinear  viscoelasticity  of  elastomers  is  commonly  described  through  a  rheological  spring- 
dashpot  models  of  the  Zener  type  (e.g.  Roland,  1989,  Johnson  et  ah,  1995,  Bergstrom  and 
Boyce,  1998,  Huber  and  Tsakmakis,  2000,  Quintavalla  and  Johnson,  2004,  Bergstrom  and 
Hilbert,  2005,  Qi  and  Boyce,  2005,  Areias  and  Matous,  2008,  Hoo  Fatt  and  Ouyang,  2008, 
Tomita  et  al.,  2008). 

In  the  present  work,  we  present  the  time  integration  scheme  for  a  newly  developed  rate- 
dependent  constitutive  model  for  polyurea  (Shim  and  Mohr,  2009b).  After  implementing  the 
model  as  a  user  material  subroutine  into  a  commercial  finite  element  software,  the  model  is  used 
to  predict  the  mechanical  response  of  thick  polyurea  layers  under  punch  loading.  Experiments 
are  performed  on  10mm  thick  polyurea  layers  for  different  punch  velocities  and  different  hemi¬ 
spherical  punch  radii.  It  is  found  that  the  model  provides  an  accurate  description  of  the  loading 
phase,  which  validates  the  assumptions  made  with  respect  to  strain-rate  and  pressure  sensitiv¬ 
ity.  However,  the  predicted  response  deviates  from  the  experimental  result  during  unloading 
which  is  discussed  in  detail. 


4.2  Punch  Experiments 

4.2.1  Specimens 

The  polyurea  specimens  used  in  the  study  are  extracted  from  a  steel  armor  plate  with  a  12.7mm 
thick  layer  of  polyurea  DragonShield-HT  Explosive  Resistant  Coating  (ERC).  Rectangular  sam¬ 
ples  of  46  x  40mm  are  cut  from  the  coated  armor  plate  using  conventional  machining.  The  coated 
polyurea  is  not  separated  from  the  steel  as  the  steel  serves  as  specimen  support  throughout  the 
punching  experiments.  However,  to  guarantee  a  uniform  layer  thickness  for  all  specimens,  the 
outer  surface  of  the  polyurea  is  machined  down  to  a  final  polyurea  thickness  of  10mm.  All 
experiments  are  performed  on  polyurea  in  its  virgin  state  (no  prior  loading)  after  a  shelf  life  of 
about  four  years. 
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4.2.2  Experimental  Procedure 


The  specimens  are  clamped  on  the  table  of  a  hydraulic  testing  machine  (Model  8080,  Instron). 
The  specimens  are  loaded  through  hemispherical  indenters  that  are  attached  to  the  moving 
actuator  of  the  upper  crosshead.  Two  hemispherical  indenters  of  different  size  are  employed: 
D  =  12.7 mm  and  D  =  44.45mm  (see  Figures  4-la  and  4-lb).  In  addition,  we  consider  two  dif¬ 
ferent  types  of  lateral  boundary  conditions:  (1)  free  in  all  lateral  directions,  and  (2)  constrained 
in  the  width  direction  (Figures  4-lc  and  4-ld).  For  the  latter  case,  the  polyurea  specimen  is 
placed  between  two  steel  blocks  which  prevents  bulging  in  the  width  direction,  but  does  not 
prevent  possible  shrinking  of  the  specimen  in  width  direction.  The  friction  at  the  interface 
between  the  indenters  and  the  polyurea  is  reduced  by  grease  and  multiple  0.1mm  thick  Teflon 
layers  (which  are  partially  torn  apart  during  the  test).  The  experiments  are  performed  under 
displacement  control  at  constant  deceleration  using  the  control  software  MAX  (Instron,  Can¬ 
ton).  Starting  with  an  initial  velocity  vq,  the  velocity-time  profile  decreases  linearly  until  the 
experiment  is  stopped  at  a  velocity  of  — t>o  (Figure  4-2).  The  initial  position  of  the  actuator 
is  chosen  such  that  the  loading  direction  is  reversed  (point  vq  —  0)  when  the  punch  depth 
reaches  7mm.  Experiments  are  performed  for  vq  =  1mm /s  and  vo  =  100 mm/s.  Throughout 
all  experiments,  the  punch  force  is  measured  using  a  50kN  load  cell.  The  punch  displacement 
is  measured  using  an  LVDT  that  is  integrated  in  the  actuator.  The  overall  stiffness  of  the  test¬ 
ing  frame  of  about  100 kN/mm  has  obtained  from  the  comparison  of  the  LVDT  readout  with 
optical  displacement  measurements  (Shim  and  Mohr,  2009a).  All  displacements  reported  in  the 
following  have  been  corrected  by  the  deformation  associated  with  the  finite  machine  stiffness. 


4.2.3  Experimental  Results 

Figure  4-3a  summarizes  the  results  for  the  experiments  with  the  small  hemispherical  indenter 
(D  =  12.7 mm).  The  measured  force-displacement  curves  are  monotonically  increasing  up  to 
the  point  of  load  reversal.  During  unloading,  the  force  reaches  zero  at  a  displacement  of  about 
3.2 mm.  Beyond  that  point,  the  punch  moves  faster  than  the  surface  of  the  creeping  polyurea 
specimen.  The  force  level  is  about  40%  higher  in  a  punch  test  at  vo  =  100 mm/s  than  that 
at  vq  —  1mm/ s  which  is  mostly  due  to  the  rate  dependent  material  behavior.  Note  that  the 
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Figure  4-1:  Photos  of  the  experiments  with  (a)  the  small  punch  and  (b)  the  large  punch. 
The  second  row  shows  the  set  up  for  (c)  free  boundaries  in  all  lateral  directions  and  (d)  for 
constrained  boundary  conditions  in  the  width  direction. 


Figure  4-2:  Applied  loading  profiles:  (a)  Applied  velocity  history;  the  velocity  axis  is  normalized 
by  the  applied  initial  velocity  of  either  wo  =  1mm/ s  or  100 m/s;  the  time  axis  is  normalized  by 
|2uo/vo|  with  u o  =  7mm;  (b)  Corresponding  applied  displacement  history. 
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Displacement  [mm]  Displacement  [mm] 

Figure  4-3:  Measured  load-displacement  curves  for  experiments  with  (a)  the  small  punch,  (b) 
the  large  punch. 

loading  conditions  are  of  quasi-static  nature  since  the  speed  of  elastic  waves  is  polyurea  (of 
the  order  of  106mm/s)  is  much  faster  than  the  loading  velocities.  It  is  interesting  to  observe 
the  convergence  of  the  force-displacement  curves  for  both  velocities  upon  unloading  which  is 
consistent  with  the  results  from  uniaxial  compression  experiments  (Shim  and  Mohr,  2009b). 
The  results  with  and  without  constraint  in  the  width  direction  are  almost  identical.  Thus, 
Figure  4-3a  presents  the  results  for  free  lateral  boundary  conditions  only. 

In  close  analogy  with  the  results  from  the  small  punch,  the  measured  force-displacement 
curves  for  the  large  hemispherical  punch  experiments  are  loading  velocity  sensitive  and  show 
a  higher  force  level  for  the  high  loading  velocity  (Figure  4-3b).  Moreover,  the  large  punch 
experiments  also  show  the  characteristic  convergence  of  the  force-displacement  curves  during 
unloading.  For  the  large  punch,  the  force  level  is  zero  at  a  punch  displacement  of  3.4mm.  The 
effect  of  the  boundary  condition  in  width  direction  becomes  apparent  when  using  the  large 
hemispherical  punch  (D  —  44.45mm);  the  force  level  with  the  constraint  (solid  lines)  is  higher 
than  that  with  the  free  lateral  boundary  conditions  (dotted  lines).  Regardless  of  the  applied 
velocity  profiles,  the  constraint  in  the  width  direction  increases  the  force  level  by  about  3 kN 
at  the  maximum  punching  depth.  Subtracting  the  force  level  with  vo  =  1mm /s  from  the  force 
level  from  vq  =  100 mm/ s  for  the  corresponding  displacement  magnitude,  one  can  find  that  the 
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Figure  4-4:  Rheological  model  of  the  rate  dependent  constitutive  model  for  polyurea. 

effect  strain  rate  during  loading  phase  results  in  the  increase  of  force  level  up  to  7kN  about  at 
the  punching  depth  of  7 mm. 

4.3  Constitutive  Model 

In  the  following,  we  present  the  algorithmic  version  of  a  recently-developed  constitutive  model 
for  polyurea.  The  reader  is  referred  to  Shim  and  Mohr  (2009b)  for  details  on  the  differen¬ 
tial  formulation  and  the  underlying  physical  arguments  for  specific  constitutive  equations. 
The  algorithm  is  implemented  as  a  user  material  subroutine  for  the  finite  element  software 
Abaqus /explicit.  The  constitutive  equations  are  based  on  a  rheological  model  of  two  Maxwell 
elements  that  act  in  parallel  (Figure  4-4).  The  first  Maxwell  element  represents  the  soft  part 
( Network  A)  of  polyurea  and  is  composed  of  a  nonlinear  viscous  damper  and  a  nonlinear  Gent 
spring.  The  hard  part  ( Network  B )  is  represented  by  another  nonlinear  viscous  damper  and 
a  Hencky  spring.  The  equations  are  cast  in  a  framework  of  finite  strains  with  multiplicative 
Kroner-Lee  decomposition  (Kroner,  1960,  and  Lee,  1969)  of  the  deformation  gradient  for  each 
Maxwell  element.  The  specific  evolution  laws  for  individual  model  components  are  given  in 
algorithmic  form  below.  The  Euler  forward  numerical  integration  method  is  employed  within 
the  subroutine  for  the  effectiveness  of  the  implementation,  because  of  numerical  challenges  as¬ 
sociated  with  Euler  backward  integration  schemes  for  Maxwell  model  with  nonlinear  rubber 
spring  and  nonlinear  dashpot  (e.g.  Areias  and  Matous,  2008). 

The  material  model  comprises  the  viscous  deformation  gradients  of  Networks  A  and  B , 
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and  FVB,  as  internal  state  variables.  We  consider  the  strain  driven  time  integration  problem, 
where  the  variables  at  time  t  —  t  +  At  are  calculated  based  on  the  solution  at  time  t.  In  other 
words,  given  the  total  deformation  gradient  F tot  (t)  and  the  internal  state  variables  (f)  and 
Fvb  ( t ),  we  evaluate  the  total  Cauchy  stress  T  (r)  along  with  the  updated  state  variables  F^  (r) 
and  Fj3  (r).  The  hydrostatic  part  of  the  Cauchy  stress  tensor  is  directly  related  to  the  change 
in  total  volume, 


J  (t)  =  det  F tot  (r) 
trT  (t)  In  J  (t) 


(4.1) 

(4.2) 


where  k  denoting  the  bulk  modulus.  The  deviatoric  part  of  the  macroscopic  stress  tensor  corre¬ 
sponds  to  the  sum  of  the  deviatoric  stresses  and  T b  acting  on  Network  A  and  Network  B, 
respectively, 


T  (r)  -  devT  (r)  +  tlT^  1 

=  T,(r)  +  TB(r)  +  ®l 


(4.3) 


For  the  evaluation  of  the  stresses  and  internal  variables  in  Networks  A  and  B,  it  is  useful  to 
define  the  isochoric  deformation  gradient 


F(t)  =  [J(t)]  1/3  Ftot  (r) 


(4.4) 


as  well  as  the  strain-like  deformation  measure 

c  (r)  =  yj tr  {FT  (r)F(r)}-3  (4.5) 

4.3.1  Response  of  Network  A 

The  viscous  gradient  is  calculated  from  the  Euler  forward  form 

F)4  (r)  =  [1  +  At  (t)]  F\  (t)  (4.6) 
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where  DVA  (t)  is  the  rate  of  viscous  deformation  at  time  t.  Subsequently,  we  have  the  isochoric 
deformation  gradient  of  elastic  deformation  in  Network  A 


Fe4(r)=F4  (t)  [FUf)]-1 


(4.7) 


Using  (r),  we  determine  the  deviatoric  Cauchy  stress  for  Network  A  based  on  Gent’s  (1996) 
free  energy  function 


T\4  (r)  - 


J  (r) 


tr{Ff  (r)FeA(r)}-3A 
Ja  J 


dev{F^(r)Ff  (t)} 


(4.8) 


with  the  material  parameters  nA  >  0  (initial  modulus)  and  JA>  0  (locking  stretch). 

The  viscous  rate  of  deformation  tensor  D'^  (r)  is  obtained  from  the  nonlinear  viscous  evo¬ 
lution  law.  For  this,  we  calculate  the  driving  Mandel  stress 


Ma  (r)  =  J  (r)  dev  {f£/  (r)  TA  (r)  F^r  (r)} 


(4.9) 


along  with  the  corresponding  equivalent  stress 


mA  (t) 


2  Ma  (t)  :  Ma  (t) 


(4.10) 


The  equivalent  rate  of  viscous  deformation  dA  (r)  is  then  given  by  the  power-law 


with  the  reference  rate  of  deformation  dref  =  1/s  and  the  material  properties  PA  >  0  (viscosity 
constant)  and  the  exponent  nA  >  0  (constant  for  rate  sensitivity).  The  flow  rule  assumes  that 
the  rate  of  viscous  deformation  tensor  is  aligned  with  the  driving  Mandel  stress  Ma 


3  dA(r) 


2  mA  (r) 


mai 
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4.3.2  Response  of  Network  B 

In  close  analogy  with  the  procedure  for  Network  A,  the  constitutive  equations  for  Network  B 
are  solved  numerically.  We  have 

F£(r)==[l  +  AtDMi)]FMi)  (4-12) 


and 

Fs  (r)  =  Fb  (t)  [Fg  (t)}-1  (4.13) 

The  stiff  elastic  response  of  Network  B  is  described  by  Hencky’s  strain  energy  function.  Thus, 
elastic  right  stretch  tensor  TJB  (r)  and  the  rotation  tensor  R|j  (r)  are  calculated  from  the  polar 
decomposition  of  the  elastic  deformation  gradient, 

F%  (r)  =  Rg  (r)  V%  (r)  with  Rf  (r)  (r)  =  1  (4.14) 


before  calculating  the  deviatoric  Cauchy  stress 

T b  (r)  -  (r)  [InUf,  (r)]  Rf  (r)  (4.15) 

with  the  shear  modulus  nB. 

Subsequently,  we  write 


Mb  (r)  =  J  (r)  dev  {Ff  (r)  TB  (r)  Ff T  (r) } 


niB  (r)  =  \/ -Mg  (r)  :  Mb  (t) 


The  rate  of  viscous  deformation  is  approximated  by 


\V3ij,bQb  \ 


<  i  —  exp 


C  (r)  (  dB  (t) 


Qb  V  dref 


-nB 


-1\  1/»B 


(4.16) 

(4.17) 


(4.18) 


where  QB  and  nB  are  material  model  parameters  that  control  the  rate  sensitivity  of  Network  B. 
As  for  Network  A,  the  rate  of  viscous  deformation  tensor  is  aligned  with  the  driving  Mandel 
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Table  4.1:  Summary  of  material  parameters  indentified  from  monotonic  loading/unloading  tests 
under  five  different  strain  rates. _ 


Rubber 

Ha 

imMPa 

Isochoric  Part 

Spring 

Ja 

10.7 

of  Network  A 

Viscous 

Pa 

4.42MPa 

Damper 

nA 

0.0646 

Linear 

Spring 

Hb 

82.3  MPa 

Isochoric  Part 
of  Network  B 

Viscous 

Qa 

0.0447 

Damper 

nA 

0.0755 

Volumetric  Part 

K 

829 MPa 

stress  Mb 

ds<tH£mMb(t)  (4J9) 

4.3.3  Model  Parameter  Identification 

The  constitutive  model  requires  the  identification  of  eight  material  parameters:  four  parameters 
Jj\ ,  PA ,  nA)  for  Network  A,  three  parameters  (  fiB,  Qa,  nA)  for  Network  B,  and  one 
parameter  (k)  describing  the  elastic  volumetric  response.  All  model  parameters  have  been  iden¬ 
tified  based  on  the  results  from  uniaxial  compression  experiments  at  five  different  strain  rates 
between  10~3/s  and  llT/s  up  to  a  true  strain  of  —1.0.  Details  on  the  parameter  identification 
procedure  are  given  in  Shim  and  Mohr  (2009b).  Table  4.1  summarizes  the  identified  material 
parameters  which  are  used  in  the  present  structural  validation  study.  The  comparison  of  the 
measured  and  predicted  stress-strain  curves  for  uniaxial  compression  is  shown  in  Figure  4-5.  It 
is  noted  that  the  present  choice  of  material  model  parameters  provides  a  good  description  of  the 
initial  monotonic  loading  response  of  polyurea,  while  the  model  predictions  are  systematically 
too  stiff  during  unloading. 
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Figure  4-5:  Comparison  of  simulation  results  and  experiments  for  continuous  loading-unloading 
cycles,  (a)  10-3/s,  (b)  10_2/s,  (c)  10_1/s,  (d)  10°/s,  (e)  101/s. 
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4.4  Numerical  Simulations  of  the  Punch  Experiments 

Finite  element  simulations  are  performed  of  all  punch  experiments.  We  make  use  of  the  sym¬ 
metry  of  the  mechanical  system  by  using  a  quarter  model.  The  polyurea  block  is  meshed  with 
eight-node  reduced  integration  solid  elements  (type  C3D8R  of  the  Abaqus  element  library) 
while  using  the  user  material  option  to  describe  the  constitutive  behavior.  The  punches  are 
modeled  using  rigid  elements.  The  meshes  comprise  eight  elements  in  thickness  direction  and 
a  small  bias  in  the  horizontal  plane  providing  a  smaller  element  size  near  the  center  than  at 
the  specimen  boundaries.  Frictionless  interface  conditions  are  assumed  between  the  punch  and 
the  polyurea.  The  motion  of  all  nodes  at  the  bottom  of  the  polyurea  (which  corresponds  to  the 
interface  with  the  steel  substrate  in  the  experiment)  is  set  to  zero.  We  applied  the  same  punch 
velocity  histories  as  in  the  experiments  and  omitted  addition  of  bulk  viscosity. 

The  deformed  meshes  at  the  point  of  maximum  penetration  (uq  =  7 mm)  are  shown  in  Figure 
4-6  for  the  simulation  without  lateral  constraint  along  the  width  direction.  The  comparison  of 
Figures  4-6a  and  4-6b  clearly  shows  that  the  bulging  effect,  i.e.  the  lateral  expansion  of  the 
polyurea  block  is  more  pronounced  for  the  large  than  for  the  small  punch.  In  the  latter  case, 
we  observe  a  small  lateral  displacement  of  about  0.5 mm.  The  bulging  has  little  effect  in  the 
load-displacement  curves  for  the  small  punch  and  we  found  nearly  the  same  force-displacement 
curve  for  free  and  constrained  lateral  boundary  conditions. 

The  predicted  force-displacement  curves  for  the  small  punch  are  shown  in  Figure  4-7a. 
The  simulations  results  are  in  good  agreement  with  the  measured  force-displacement  curves 
during  the  loading  phase.  A  similar  conclusion  may  be  drawn  from  comparing  the  simulations 
and  experiments  for  the  large  punch  (Figure  4- 7b  and  4-7c).  The  comparison  shows  that  the 
numerical  model  can  predict  (a)  the  effect  of  the  punch  size,  (b)  the  effect  of  the  lateral  boundary 
constraint,  and  (c)  the  effect  of  loading  velocity  on  the  force-displacement  curve  during  loading. 
Figure  4-8  shows  the  histories  of  the  strain-like  variable,  (,  and  the  viscous  rates,  d/ \  and 
dB,  for  the  simulation  of  a  small  punch  with  i.’o  =  100 mm/s  (the  corresponding  locations  are 
highlighted  in  Figure  4-6a).  The  profiles  elucidate  the  strong  variations  in  local  strain  rate 
during  the  punch  experiments.  It  is  interesting  to  see  that  the  local  viscous  rate  can  be  as  high 
as  30/s  which  is  close  to  the  strain  rate  of  the  fastest  calibration  experiment  (Figure  4-5e). 
However,  unlike  for  the  calibration  experiments,  the  strain  rates  are  non-constant  throughout 
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Figure  4-6:  The  contour  plots  of  the  logarithmic  strain  in  thickness-direction  from  simulations 
with  vo  —  100 m/s  at  an  indentation  depth  of  7mm:  (a)  small  hemispherical  indenter,  and  (b) 
large  hemispherical  indenter.  The  detail  in  Figure  4-6a  shows  the  locations  for  which  the  strain 
rates  d a  and  dg  are  plotted  in  Figure  4-8. 
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a  punch  tests  which  is  seen  as  an  important  validation  of  the  model  assumptions  with  respect 
to  the  effect  of  strain  rate.  At  a  punching  depth  of  7 mm,  the  simulation  results  overestimate 
the  force  level  by  about  0.5 kN  in  the  case  of  the  small  punch  and  2.5 kN  in  the  case  of  the 
large  punch.  This  overestimation  (which  corresponds  to  less  than  10%  of  the  current  force 
level)  is  consistent  with  the  model  calibration  results.  Recall  from  Figure  4-5  that  the  model 
systematically  overestimates  the  stress  level  at  very  large  compressive  strains. 

The  model  predictions  during  unloading  deviate  from  the  experimentally-measured  force- 
displacement  responses.  The  simulations  predict  the  characteristic  convergence  of  the  load- 
displacement  curves,  but  the  predicted  response  during  unloading  is  too  stiff  for  both  velocities 
and  both  punch  sizes.  Consequently,  the  instantaneous  residual  displacement  at  zero  force 
is  overestimated  in  the  simulations.  This  deviation  is  again  consistent  with  the  calibration 
experiments,  indicating  an  inherent  shortcoming  of  the  model  formulation  for  unloading. 

4.5  Discussion  on  Unloading  Behavior  of  Polyurea 

The  loading  response  of  the  model  agrees  well  with  the  test  results  except  for  an  overestimation 
of  the  force  at  large  punch  displacements.  The  simulation  results  indicate  that  the  material  is 
subject  to  strains  of  up  to  —2.0  as  the  punching  depth  reaches  7 mm  (see  Figure  4-6).  Since 
this  is  twice  as  high  as  in  the  experiments  for  the  material  model  parameter  identification,  we 
performed  additional  material  tests  at  a  strain  rate  of  10~2/s.  Figure  4-9a  shows  the  resulting 
stress-strain  curve  from  three  different  experiments.  Each  virgin-state  specimen  is  subject  to  a 
single  loading- unloading  cycle  up  to  a  maximum  strain  of  -0.5,  -1.0  and  -1.5,  respectively.  The 
corresponding  simulations  (Figure  4-9b)  reveal  that  the  numerically-predicted  locking  behavior 
is  more  severe  than  that  observed  in  the  experiments. 

The  unloading  path  in  the  stress-strain  curves  comprises  two  characteristic  regimes:  a  stiff 
part  at  the  beginning  of  unloading  followed  by  a  soft  part  as  the  stress  approaches  zero.  Both 
the  experiments  and  the  numerical  simulations  exhibit  this  feature.  The  numerical  model  does 
not  show  the  smooth  transition  between  these  two  regimes,  but  this  is  seen  as  an  acceptable 
engineering  approximation  of  the  physical  behavior.  However,  the  apparent  main  deficiency  of 
the  current  model  is  its  inability  to  capture  the  increase  of  the  hysteresis  loop  width  as  the  strain 
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Figure  4-7:  Comparison  of  simulations  and  experiments.  Force-displacement  curves  for  (a) 
small  punch  with  free  lateral  boundaries,  (b)  large  punch  with  the  free  lateral  boundaries,  (c) 
large  punch  with  constraint  in  the  width  direction. 
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Figure  4-8:  Results  from  the  small  punch  simulation  with  vo  =  100 mm/s  and  free  lateral 
boundary  conditions  at  three  different  locations  (labeled  by  A,  B  and  C  in  Figure  4-6(a): 
Histories  of  (a)  the  strain-like  deformation  measure,  £,  (b)  the  viscous  strain  rate  of  Network  A, 
d,A,  (c)  the  viscous  strain  rate  of  Network  B,  dg-  Note  that  the  loading  direction  is  reversed 
at  t  —  0.14s. 
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Figure  4-9:  Stress-strain  curves  for  a  single  loading-unloading  cycle  at  e  =  10  2 / s  as  obtained 
from  (a)  experiments  and  (b)  simulations, 

increases.  As  the  applied  maximum  strain  increases,  the  simulations  predict  almost  constant 
stress  drops  (marked  by  a1,  b'  and  d  in  Figure  4-9b)  while  the  experiments  show  a  substantial 
increase  in  the  magnitude  of  the  stress  drops  as  a  function  of  strain  (marked  by  a,  b  and  c  in 
Figure  4-9a). 

Recall  that  Network  A  in  the  current  constitutive  model  is  mainly  responsible  for  the 
rubbery  behavior  while  Network  B  describes  the  high  initial  stiffness  and  the  time-dependent 
hysteresis.  The  contributions  of  Networks  A  and  B  are  in  opposite  direction  as  far  as  the 
hysteresis  width  is  concerned.  To  shed  more  light  on  this  particular  feature,  Figure  4-10  shows 
a  direct  comparison  of  the  model  response  and  the  experiments.  In  addition,  we  plotted  the 
individual  contributions  of  Network  A  and  Network  B  to  the  stress-strain  curve  (dashed  curves 
in  Figure  4-10).  Note  that  the  magnitude  of  the  stress  contribution  of  Network  B  depends  on 
the  compressive  strain  only,  while  its  sign  changes  from  compression  to  tension  upon  unloading. 
The  contribution  of  Network  A  on  the  other  hand  is  a  compressive  stress  irrespective  of  the 
loading  direction.  Thus,  as  the  compressive  strain  increases,  Network  B  makes  the  hysteresis 
wider  while  the  opposite  holds  true  for  Network  A.  In  order  to  replicate  the  experimental 
observation  of  a  hysteresis  width  increase  as  a  function  of  the  maximum  compressive  strain 
(Figure  4-9a),  the  contribution  of  Network  B  should  be  dominant.  Moreover,  the  contribution 
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of  Network  B  would  need  to  increase  as  the  compressive  strain  increases.  However,  due  to  the 
specific  choice  of  the  viscous  evolution  law  for  Network  B,  the  stress  contribution  of  Netivork  B 
is  more  or  less  constant  once  the  compressive  strain  has  exceeded  —0.1  (see  the  stress  plateau 
in  the  response  curves  for  Network  B  in  Figure  4-10).  An  attempt  was  made  to  change 
the  calibration  of  the  response  of  Network  B  using  the  present  modeling  framework,  but  the 
subsequent  simulation  results  were  no  longer  satisfactory  for  the  phase  of  loading. 

A  constitutive  model  with  a  different  rheological  composition  needs  to  be  used  to  improve 
the  predictions  for  unloading.  A  previous  study  (Shim  and  Mohr,  2009b)  has  shown  that 
models  of  the  Zener-type  (spring  in  parallel  with  a  single  Maxwell  element)  cannot  describe 
the  large  deformation  behavior  of  polyurea  over  a  wide  range  of  strain  rates.  The  present 
results  indicate  that  the  assumption  of  two  Maxwell  elements  in  parallel  provides  an  accurate 
description  for  monotonic  loading  only.  It  can  describe  the  two  characteristic  stiff  and  soft 
regimes  during  unloading  over  a  wide  range  of  strain  rates,  but  the  model  predictions  are  only 
in  poor  quantitative  agreement  during  unloading. 

Throughout  our  model  development,  we  focused  on  a  single  loading-unloading  cycle  on 
virgin-state  specimens.  However,  even  though  the  loading  of  the  material  in  its  virgin  state 
appears  to  be  the  most  important  with  respect  to  the  real-life  applications  of  polyurea,  the 
material  response  to  cyclic  loading  may  be  instructive  as  far  as  the  choice  of  the  rheological 
model  is  concerned.  Figure  4-11  shows  the  stress-strain  curve  for  polyurea  for  five  consecutive 
loading-unloading  cycles.  After  each  loading-unloading  cycle  at  a  constant  true  strain  rate  of 
10~"2/s,  we  let  the  specimen  creep  at  zero  stress  for  about  300s  before  applying  the  subsequent 
loading-unloading  cycle.  The  significant  difference  between  the  first  and  subsequent  loading 
cycles  illustrates  the  Mullins  effect  for  polyurea.  The  stress-strain  curve  for  the  first  loading  is 
characterized  by  the  high  initial  stiffness  and  the  high  peak  stress  marked  as  A.  After  the  first 
loading-unloading  cycle,  however,  the  shape  of  the  loading  path  changes  noticeably  while  nearly 
the  same  unloading  path  is  observed.  Additional  loading-unloading  cycles  create  very  little 
changes  in  both  loading  and  unloading.  This  observation  suggests  that  the  constitutive  model 
should  comprise  an  internal  variable  that  reflects  the  amount  of  “microstructural  damage” 
associated  with  the  Mullins  effect.  The  Mullins  effect  has  been  investigated  by  many  research 
groups  using  either  damage-based  constitutive  models  (e.g.,  Simo,  1987,  Govindjee  and  Simo, 
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Figure  4-10:  Comparison  of  simulation  results  and  experiments  for  continuous  load¬ 
ing/unloading  cycles  at  e  =  10 ~2/s.  (a)  £/  =  —0.5,  (b)  £y  =  —1.0,  (c)  £y  =  —1.5. 
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Tests  at  10'2/sec 


Figure  4-11:  Illustration  of  the  Mullins  effect.  Five  compression  loading  and  unlading  cycles 
are  performed  at  the  constant  strain  rate  of  10-2/s  up  to  the  maximum  strain  of  —1.0.  The 
stress  is  zero  between  subsequent  cycles  for  about  five  minutes. 

1991,  1992,  Lion,  1996,  1997,  Miehe  and  Keck,  2000)  or  the  concept  of  hard/soft  domain 
reorganization  (e.g.  Johnson  and  Beatty,  1993a,  1993b,  Beatty  and  Krishnaswamy,  2000,  Qi 
and  Boyce,  2004,  2005).  Here,  the  explicit  account  of  the  Mullins  effect  is  deferred  to  future 
work  since  further  experimental  data  is  needed  to  analyze  the  effect  of  loading  velocity  over  a 
wide  range  of  strain  rates. 

4.6  Conclusion 

Punch  experiments  have  been  performed  on  10mm  thick  polyurea  samples  at  punch  velocities 
of  up  to  100 mm/s  which  resulted  in  maximum  local  viscous  strain  rates  of  up  to  30/s  inside  the 
polyurea  layers.  A  newly-developed  rate-dependent  constitutive  model  for  polyurea  has  been 
implemented  into  a  finite  element  program  and  used  to  predict  the  experimentally-measured 
force-displacement  curves  for  different  punch  sizes  and  velocities.  Furthermore,  the  effect  of 
lateral  constraints  has  been  investigated.  The  model  provides  accurate  predictions  of  the  load¬ 
ing  response  for  all  six  test  configurations  which  is  interpreted  as  a  partial  validation  of  the 
assumptions  made  with  respect  to  strain-rate  and  pressure  sensitivity.  During  unloading,  the 
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model  exhibits  the  characteristic  stiff  and  soft  responses  of  polyurea;  however,  it  systematically 
overestimates  the  depth  of  the  quasi-instantaneous  residual  punch  imprint  after  unloading. 
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Chapter  5 


Conclusion  and  Suggestions 


5.1  Summary  of  Main  Results 

The  main  results  of  this  dissertation  can  be  summarized  twofold:  experiments  and  constitutive 
modeling.  First,  the  effect  of  strain  rate  on  polyurea  are  experimentally  quantified  for  a  wide 
range  of  strain  rates,  i.e.  10~3/s  ~  4000/s.  The  modified  SHPB  system  of  Zhao  and  Gary 
(1997)  has  been  used  to  perform  compression  tests  on  polyurea  at  low,  intermediate  and  high 
strain  rates  (10/s  ~  1000/s).  It  is  composed  of  nylon  input  and  output  bars,  while  the  striker 
bar  is  substituted  by  a  hydraulic  actuator.  Using  the  deconvolution  technique  by  Bussac  et  al. 
(2002),  the  time  limitation  of  conventional  SHPB  systems  may  be  overcome,  thereby  enabling 
the  use  of  the  modified  SHPB  system  for  low  and  intermediate  strain  rate  experiments  of 
long  duration.  The  experiments  confirm  the  known  strain  rate  sensitivity  of  polyurea.  The 
measured  stress  levels  correspond  well  to  earlier  results  which  have  been  obtained  from  tests 
on  conventional  SHPB  systems  with  very  long  bars. 

Second,  a  new  constitutive  model  is  proposed  to  capture  the  rate-sensitivity  for  a  wide  range 
of  strain  rates  and  the  characteristic  high  stiffness  behavior  of  virgin  polyurea  at  small  strains. 
Various  types  of  loading  conditions  are  tested  for  the  purpose  of  validation,  and  their  results 
are  well  compared  with  the  model  predictions  in  loading  phases.  The  large  strain  compression 
response  of  polyurea  is  investigated  for  strain  rates  ranging  10~3/s  from  to  101/s.  Continuous 
and  stair-like  loading  and  unloading  experiments  are  performed  in  addition  to  simple  relaxation 
tests.  The  experimental  results  reveal  that  the  so-called  equilibrium  path  concept  breaks  down 
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in  the  case  of  polyurea.  In  contradiction  with  the  definition  of  the  equilibrium  path,  the  iden¬ 
tified  average  stress-strain  curves  in  continuous  loading/unloading  experiments  depend  on  the 
rate  of  loading.  Thus,  as  an  alternative  to  the  equilibrium  path  based  rheological  models  of  the 
Zener  type,  a  new  constitutive  model  is  formulated  assuming  two  parallel  Maxwell  elements. 
The  finite  strain  constitutive  equations  are  outlined  in  detail,  and  it  has  been  implemented 
into  a  finite  element  program.  Subsequently,  the  eight  material  model  constants  are  calibrated 
to  describe  the  mechanical  behavior  of  polyurea.  The  model  predictions  are  in  good  agree¬ 
ment  with  the  experimental  results.  The  characteristic  convergence  of  the  unloading  paths  for 
different  strain  rates  is  successfully  captured  by  the  model.  Furthermore,  the  rate  dependent 
viscous  dissipation  throughout  a  loading/unloading  cycle  of  uniaxial  conditions  is  described 
with  reasonable  accuracy.  For  the  validation  of  the  proposed  model,  punch  experiments  have 
been  performed  on  10mm  thick  polyurea  samples  at  punch  velocities  of  up  to  100 mm/s  which 
resulted  in  maximum  local  viscous  strain  rates  of  up  to  30/s  inside  the  polyurea  layers.  The 
proposed  constitutive  model  is  used  to  predict  the  experimentally-measured  force-displacement 
curves  for  different  punch  sizes  and  velocities.  Furthermore,  the  effect  of  lateral  constraints  has 
been  investigated.  The  model  provides  accurate  predictions  of  the  loading  response  for  all  six 
test  configurations  which  is  interpreted  as  a  partial  validation  of  the  assumptions  made  with 
respect  to  strain-rate  and  pressure  sensitivity. 

5.2  Suggestions  for  Future  Studies 

The  following  topics  can  be  suggested  for  future  studies: 

•  Constant  strain  rates  from  the  hydraulic  SHPB:  Although  the  intrinsic  time  limitation 
of  SHPB  systems  could  be  overcome,  this  study  also  shows  that  it  is  still  not  possible 
to  perform  experiments  at  reasonably  constant  strain  rates  with  this  technique.  This  is 
due  to  the  finite  length  (i.e.  3m  in  this  study)  of  the  input  and  output  bars  which  causes 
a  periodic  change  in  loading  velocity.  It  is  shown  that  intermediate  strain  rate  SHPB 
experiments  require  either  very  long  bars  (>  20m)  or  very  short  bars  (<  0.5m)  in  order 
to  achieve  an  approximately  constant  strain  rate  throughout  the  entire  experiment. 

•  Improvement  in  the  fixed  boundary  conditions  in  the  hydraulic  SHPB :  In  order  to  prevent 
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the  failure  of  the  nylon  bars  under  excessive  loads,  a  fixed  end  support  system  in  the 
hydraulic  SHPB  of  the  current  study  is  designed  such  that  the  bars  are  released  before 
elastic  buckling  occurs.  However,  with  the  current  design  of  the  fixed  end,  significant 
force  drops  are  found  during  tests  due  to  the  premature  partial  failure  of  the  fixed  end 
support  of  the  output  bar.  A  new  design  of  the  fixed  end  support  can  be  investigated  to 
prevent  the  elastic  buckling  of  the  bars  and  to  resist  the  specimen  strength  at  the  same 
time. 

•  Improvement  in  the  unloading  prediction  of  the  constitutive  model:  From  various  tests 
with  different  loading  conditions,  the  model  exhibits  the  characteristic  stiff  and  soft  re¬ 
sponses  of  polyurea  during  unloading.  However,  the  predicted  unloading  behavior  is  too 
stiff,  and  the  model  systematically  overestimates  the  depth  of  the  quasi-instantaneous 
residual  strain/displacement  after  unloading.  A  new  rheological  model  can  be  investi¬ 
gated  to  improve  the  prediction  of  unloading  behavior  and  to  capture  the  Mullins  effect 
for  polyurea. 
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Appendix  A 


Identification  of  the  Wave 
Propagation  Coefficient  for 
Viscoelastic  Bars 


The  complex- valued  propagation  coefficient  £  (oj)  is  a  function  of  both  the  geometric  and  ma¬ 
terial  properties  of  the  bars.  If  the  complex  modulus  of  the  viscoelastic  bar  material  is  known, 
£  (u>)  may  be  calculated  from  solving  Pochhammer-Chree’s  frequency  equation  (Zhao  and  Gary, 
1995).  As  an  alternative,  we  make  use  of  an  experimental  method  that  considers  both  geometric 
dispersion  and  viscous  attenuation  within  the  framework  of  the  1-D  wave  theory  (e.g.  Bacon, 
1998,  Lundberg  and  Blanc,  1988)  to  determine  £  ( u ).  Recall  the  solution  of  the  one-dimensional 
wave  equation  for  viscoelastic  bars, 

e  (x,  u)  =  eR  (u)  e-^x  +  eL  (w)  (A.l) 


where  £r  (x,ui)  and  £r  (x,cu)  are  the  rightward  and  the  leftward  traveling  strain  waves.  Here, 
(,  (to)  is  the  propagation  coefficient  of  the  bars  defined  by 


£  (w)  =  k  (uj)  +  ice  (w) 


LO 

c{u) 


+  ia  (w) 


(A.2) 
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where  k  (w)  is  the  wave  number,  c  (uj)  is  the  frequency-dependant  longitudinal  wave  propagation 
speed,  and  a  ( ui )  <  0  represents  the  attenuation  coefficient.  Note  that  c  (c o)  and  a  (w)  are  even 
functions  of  u  while  k  ( u> )  is  an  odd  function  in  the  frequency  space. 

After  performing  an  impact  test  on  a  single  bar,  we  measure  the  incident  and  reflected 
strain  histories  (i.e.  e,;„c(x  —  0 ,t)  and  eref  (x  =  0 ,t))  at  x  =  0  within  in  the  bar.  In  order  to 
determine  the  transfer  function  of  the  bar,  we  use  the  two  facts  from  the  impact  test: 

•  The  rightward  and  the  leftward  travelling  strain  waves  are  equal  to  the  measured  incident 
and  reflected  strain  histories,  respectively: 

eR  (u)  =  £inc  (. x  =  0,  u)  (A.3) 


£L  {(Jj)  =  iref  {x  =  0,w) 


(A.4) 


•  The  normal  force  is  zero  at  the  free  end  of  the  bar  (i.e.  the  opposite  end  to  the  impacting 
end  of  the  bar): 

eR  (cj)  e-^d  +  el  H  e^d  =  0  (A.5) 

where  x  =  d  is  the  distance  between  the  strain  gage  location  and  the  free  end  of  the  bar. 
Now,  the  transfer  function  of  the  bar  H  (ui)  is  experimentally  determined  from 

H  M  =  _^/  (*  =  0,<A)  =  e-2 i^)d  (A.6) 

£inc  {x  =  0,w) 

Prom  Eq.  (A.6),  two  components  of  the  propagation  coefficient  can  be  identified  using  the 
relations 


«M=  arg[f,(“)l 

2  a 

(A.7) 

(A.8) 
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Appendix  B 


List  of  Papers  with  Reference  to 
Respective  Chapters 


The  following  papers  are  published  or  submitted  for  publication  with  reference  to  the  respective 
chapters. 

Chapter  2 

Shim,  J.  and  Mohr,  M.  (2009)  "Using  split  Hopkinson  pressure  bars  to  perform  large  strain 
compression  tests  on  polyurea  at  low,  intermediate  and  high  strain  rates,"  International  Journal 
of  Impact  Engineering ,  36:1116-1127. 

Chapter  3 

Shim,  J.  and  Mohr,  M.  (2009)  "Finite  strain  constitutive  model  of  polyurea  for  a  wide  range 
of  strain  rates,"  (submitted  for  publication,  September  2009). 

Chapter  4 

Shim,  J.  and  Mohr,  M.  (2009)  "Punch  indentation  of  polyurea  at  different  loading  velocities: 
Experiments  and  numerical  simulations,"  (submitted  for  publication,  October  2009). 
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